NEAR-FIELD  PRESSURE  DISTRIBUTIONS  TO  ENHANCE  SOUND 
TRANSMISSION  INTO  MULTI-LAYER  MATERIALS 


A  Dissertation 
Submitted  to  the  Faculty 
of 

Purdue  University 
by 

Andrew  M.  Jessop 


In  Partial  Fulfillment  of  the 
Requirements  for  the  Degree 
of 

Doctor  of  Philosophy 


December  2013 
Purdue  University 


West  Lafayette,  Indiana 


Report  Documentation  Page 

Form  Approved 

OMB  No.  0704-0188 

Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 

VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 

1 .  REPORT  DATE  2.  REPORT  TYPE 

10  DEC  2013  Final 

3.  DATES  COVERED 

1  Aug  2010  -  31  Dec  2016 

4.  TITLE  AND  SUBTITLE 

Near-field  Pressure  Distributions  to  Enhance  Sounds  Transmission  into 
Multi-layer  Materials  /  PhD  Dissertation 

5a.  CONTRACT  NUMBER 

N/A 

5b.  GRANT  NUMBER 

N00014-10-1-0958 

5c.  PROGRAM  ELEMENT  NUMBER 

N/A 

6.  AUTHOR(S) 

Andrew  M.  Jessop  Mechanical  Engineering  Purdue  University 

5d.  PROJECT  NUMBER 

N/A 

5e.  TASK  NUMBER 

N/A 

5f.  WORK  UNIT  NUMBER 

N/A 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Ray  W.  Herrick  Laboratories  Purdue  University  West  Lafayette,  IN 

47907 

8.  PERFORMING  ORGANIZATION  REPORT 
NUMBER 

64 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

Office  of  Naval  Research  (ONR)  875  North  Randolph  Street  -  Suite  1425 
Code  03R  Arlington,  VA  22203-1995 

10.  SPONSOR/MONITOR'S  ACRONYM! S) 

ONR 

11.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 

N/A 

12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release,  distribution  unlimited 

13.  SUPPLEMENTARY  NOTES 

N/A,  The  original  document  contains  color  images. 

14.  ABSTRACT 

The  large  impedance  difference  between  air  and  most  solids  prevents  significant  energy  transfer  from 
incident  acoustic  waves  across  the  air-material  interface.  Refraction  also  plays  a  role  in  preventing  acoustic 
transmission,  as  the  wave  speed  difference  between  air  and  solid  materials  results  in  an  increase  of  the 
resulting  propagation  angles,  creating  near-field  pressure  distributions  in  the  solid  material.  By  utilizing 
evanescent  pressure  distributions,  which  decay  normal  to  the  usual  direction  of  propagation  and  are 
represented  as  plane  waves  propagating  with  complex  angles,  energy  propagation  through  the  interface 
can  be  increased  in  the  subsonic  region  of  wave  propagation:  i.e.,  where  waves  typically  do  not  propagate 
into  a  material  with  any  effectiveness.  By  using  an  array  of  sources,  it  is  possible  to  produce  evanescent 
pressure  distributions  in  the  solid.  The  way  in  which  the  characteristics  of  this  array  of  sources  affect  the 
efficiency  of  the  generation  of  evanescent  pressure  distributions  are  explored.  Because  high  impedance 
materials  can  be  paired  to  a  lower  impedance  materials  of  interest  to  impede  acoustical  energy 
transmission,  the  wave  propagation  through  multi-layer  materials  must  be  considered  to  give  an  full 
accounting  of  power  transmission  into  structures.  A  model  for  wave  propagation  in  multi-layer  systems  of 
solids  and  fluids  was  developed  using  wave  potentials  in  each  layer,  allowing  for  coupling  between  material 
types  and  calculation  of  inter-layer  states.  By  using  wavenumber-frequency  analysis,  it  is  possible  to  target 
specific  components  in  the  multi-layer  system  and  understand  the  particular  wavetypes  that  cause  energy 
propagation  into  the  system. 


15.  SUBJECT  TERMS 

16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION 

OF  ABSTRACT 

SAR 

18.  NUMBER 
OF  PAGES 

178 

19a.  NAME  OF 
RESPONSIBLE  PERSON 

a.  REPORT 

unclassified 

b.  ABSTRACT 

unclassified 

c.  THIS  PAGE 

unclassified 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std  Z39-18 


11 


To  my  grandparents,  George  and  Mary  Jessop  and  Connie  and  Marty  Mclhorn. 
When  they  stressed  the  value  of  education,  I  don’t  think  they  ever  thought  I  would 

run  with  it  quite  like  this. 


Ill 


ACKNOWLEDGMENTS 

Dr.  J.  Stuart  Bolton  served  as  my  advisor  for  this  project.  His  input,  edits,  and 
support  were  invaluable  throughout  my  time  at  Purdue. 

The  support  system,  both  intellectually  and  socially,  provided  by  my  labmates 
over  the  years  was  a  tremendous  asset.  Tyler  Dare,  Sarah  McGuire,  Andrew  Marshall, 
Frank  Eberhardt,  Ryan  Schultz,  Chide  Giacomoni,  Wong-Hong  Choi,  YangFan  Liu, 
Mike  Hayward,  Nick  Sakamoto,  Andy  McMullen,  and  Brandon  Sobecki  (and  anyone 
else  I  might  have  left  out)  all  contributed  in  some  way  to  this  document,  whether 
they  wanted  to  or  not. 

Jelcna  Paripovic,  Jake  Miller,  Chris  Watson,  and  Daniel  Woods  worked  on  the 
MLIRI  project  along  with  me.  Best  of  luck  to  them  in  completing  their  studies.  They, 
along  with  Drs.  Davies,  Rhoads,  and  Adams  provided  spectacular  feedback  through 
our  weekly  meetings  and  progress  reports. 

The  support  staff  in  the  Herrick  Laboratories  shop  and  front  office  were  a  tremen¬ 
dous  help  as  well.  In  the  shop,  Bob  Brown,  Frank  Lee,  Ron  Evans,  Fritz  Peacock, 
and  Gilbert  Gordon  assisted  in  setting  up  tests  and  maintaining  and  evaluating  equip¬ 


ment. 


IV 


TABLE  OF  CONTENTS 

Page 

LIST  OF  TABLES .  vi 

LIST  OF  FIGURES .  vii 

ABSTRACT  .  xiii 

CHAPTER  1.  INTRODUCTION .  1 

1.1  Motivation .  2 

1.2  Contribution  of  Research .  3 

1.3  Outline .  4 

CHAPTER  2.  PREVIOUS  RESEARCH .  5 

2.1  Near-Field  Effects  on  Sound  Transmission  .  5 

2.2  Wave  Propagation  in  Multi-Layer  Systems .  7 

2.3  Generation  of  Arbitrary  Pressure  Distributions .  10 

2.4  Conclusions .  12 

CHAPTER  3.  NEAR-FIELD  EFFECTS  ON  SOLID-FLUID  ACOUSTICAL 

TRANSPARENCY .  14 

3.1  Evanescent  Plane  Waves .  15 

3.1.1  Fundamentals  of  Wave  Propagation .  15 

3.1.2  Generation  of  Evanescent  Waves  Using  Complex  Angles  ...  17 

3.1.3  Pressure  and  Velocity  Phase  Differences  in  Evanescent  Waves  20 

3.2  Evanescent  Wave  Transmission  Across  An  Interface .  21 

3.2.1  Theory  and  Principles .  21 

3.2.2  Fluid-Fluid  Interface .  27 

3.2.3  Fluid-Solid  Interface .  31 

3.3  The  Monopole  as  an  Evanescent  Wave  Source .  36 

3.3.1  Theory .  37 

3.4  Experimental  Measurement  of  Monopole  Evanescent  Components  .  39 

3.4.1  Test  Setup  and  Methodology .  40 

3.4.2  Results .  44 

3.5  Conclusions .  48 

CHAPTER  4.  GENERATION  OF  NEAR-FIELD  PRESSURE  DISTRIBU¬ 
TIONS  .  51 

4.1  Properties  of  Waves  and  Sources .  52 

4.1.1  Properties  of  Plane  Waves .  52 

4.1.2  Properties  of  the  Acoustic  Monopole .  53 


V 


Page 

4.1.3  Properties  of  the  Acoustic  Dipole .  54 

4.2  Least-Squares  Approximation  of  Pressure  Distributions  using  Multiple 

Sources .  55 

4.3  Input  Parameter  Effects  on  Plane  Wave  Fit  Using  Multi-Source  Arrays  59 

4.3.1  Incident  Wave  Parameters .  60 

4.3.2  Geometry  Parameters  .  64 

4.3.3  Source  Length  and  Density .  65 

4.3.4  Modeling  Evanescent  Waves .  67 

4.3.5  Dipole  Multi-Source  Arrays .  68 

4.4  Conclusions .  71 

CHAPTER  5.  WAVE  PROPAGATION  IN  MULTI-LAYER  MATERIALS  .  73 

5.1  Multi-Layer  Wave  Propagation  Model .  73 

5.1.1  Wave  Potential  Model  Formulation .  75 

5.1.2  Model  Verification  .  86 

5.2  Wavenumber-Frequency  Decomposition  of  Vibration .  88 

5.3  Modeling  of  Wave  Propagation  .  90 

5.3.1  The  Effect  of  Evanescent  Incident  Waves .  98 

5.4  Conclusions .  103 

CHAPTER  6.  EXPERIMENTAL  MEASUREMENTS  OF  WAVE  PROPAGA¬ 
TION  .  105 

6.1  Test  Methodology .  106 

6.2  Results .  Ill 

6.2.1  Polycarbonate  Tests  .  Ill 

6.2.2  Surrogate  Tests .  114 

6.3  Conclusions .  116 

CHAPTER  7.  CONCLUSIONS .  119 

7.1  Future  Work .  122 

LIST  OF  REFERENCES .  124 

APPENDICES 

APPENDIX  A:  MULTI-LAYER  PROPAGATION  CODE  .  130 

APPENDIX  B:  MULTI-LAYER  STATE  CALCULATION  CODE  ....  147 

APPENDIX  C:  MULTI-SOURCE  CODE .  158 

VITA  .  162 


VI 


LIST  OF  TABLES 

Table  Page 

3.1  Density  and  wave  speed  of  materials  of  interest  for  impedance  calculation.  23 

4.1  Table  of  variables  applicable  to  multi-source  arrays  and  the  desired  pres¬ 
sure  distribution .  59 

5.1  Table  of  material  properties  used  in  multi-layer  wave  propagation  model.  93 

6.1  Equipment  used  in  wave  propagation  testing .  Ill 

6.2  Parameters  for  sample  testing .  Ill 


LIST  OF  FIGURES 


Figure  Page 

3.1  Plane  wave  incidence  on  a  plane  surface.  An  incoming  wave  with  k  = 
oj/c  =  27t/A  at  an  angle  6  will  have  a  trace  wavenumber  on  the  surface  of 
kx  =  ksm(O) .  16 


3.2  The  effect  of  complex  wavenumbers  on  waveform  shape.  The  undamped 

wave  is  shown  in  the  left  plot.  A  1%  damping  (center  plot)  will  cause  the 
amplitude  to  decrease  over  distance,  while  a  1%  negative  damping  (right 
plot)  will  cause  the  amplitude  to  increase  over  distance .  17 

3.3  The  effects  of  a  complex  angle  of  propagation  (left)  compared  to  complex 

wavenumber  (right).  The  left  pressure  distribution  is  a  real  wavenumber 
with  a  complex  angle,  causing  a  decay  normal  to  the  direction  of  propa¬ 
gation.  The  right  distribution  shows  a  complex  wavenumber  propagating 
along  a  real  angle;  the  decay  is  along  the  direction  of  propagation.  ...  19 

3.4  The  resulting  trace  pressure  generated  by  incident  plane  and  evanescent 
waves  upon  a  surface.  The  upper  plots  show  propagation  over  two-dimensional 


space  above  the  incident  plane,  while  the  bottom  plots  show  the  trace  pres¬ 
sure  distribution  at  the  z  =  0  surface.  While  the  trace  pressure  of  the  30 
degree  incident  plane  wave  shown  on  the  left  does  not  decay  over  distance, 
the  resulting  trace  pressure  of  the  30  +  .1  j  degree  incidence  plot  on  the 
right  decays  along  the  plane  surface .  20 

3.5  Complex  angle  refraction  of  the  transmitted  wave  6\  =  9r  —  jdi  across  an 

interface  with  ci/c0  =  10.  The  left  plot  is  5R(arcsin(ci/c0  sin#0)),  while  the 
right  plot  is  A(ci  /co  arcsin(sin  0O)) .  26 

3.6  Pressure  transmission  coefficient  T  for  plane  wave  propagation  into  a  fluid.  28 

3.7  Pressure  transmission  coefficient  T  for  fluid-fluid  system .  29 

3.8  Intensity  transmission  coefficient  TC  variation  for  a  fluid-fluid  system  with 

a  complex  angle  of  incidence .  31 

3.9  Transmission  coefficients,  201og(Xz,)  and  201og(Ts),  for  a  system  with 

wave  speed  ratio  ci/c0  =  10  and  density  ratio  pi/po  =  1000 .  33 


Vlll 


Figure  Page 

3.10  Net  normal  intensity  above  and  below  the  fluid-solid  interface,  with  the 

incident  and  reflected  waves  contributing  to  the  intensity  above  the  in¬ 
terface  and  the  longitudinal  and  shear  transmitted  waves  contributing  to 
energy  below  the  interface .  34 

3.11  Intensity  transmission  coefficients  at  1000  Hertz  across  the  interface  for 
the  compressional  (left)  and  shear  (right)  waves  due  to  a  complex  angle  of 

the  incident  wave .  35 

3.12  Intensity  transmission  coefficients  for  a  depth  of  1  meter  into  the  material 

for  the  compressional  and  shear  waves  due  to  complex  angles .  36 

3.13  Monopole  pressure  distribution  in  the  plane  containing  the  monopole, 

showing  the  real  component  of  the  pressure  for  a  1000  Hz  monopole  in 
air  (a  wavenumber  of  18.3  radians  per  meter),  in  two-dimensional  (left) 
and  in  one- dimensional  space  (right)  along  the  radius.  The  maximum 
pressure  peak  at  r  =  0  has  been  cut  off  at  10  Pa  to  show  detail .  37 

3.14  Trace  pressure  distributions  (real  component  of  pressure)  of  a  monopole  at 

a  Zo  distance  of  10  and  200  centimeters  for  a  wavenumber  of  18.3  radians 
per  meter.  The  shorter  standoff  distance  displays  similar  properties  to  an 
evanescent  wave .  38 

3.15  Wavenumber  spectra  of  the  trace  pressure  distribution  from  a  monopole 

radiating  at  100  Hz  for  standoff  distances  of  10  and  200  centimeters.  The 
peak  in  the  plot  is  at  ko]  energy  above  this  peak  is  subsonic,  with  the 
closer  standoff  having  more  subsonic  energy  than  the  further  standoff.  .  40 

3.16  Ratio  Q  of  evanescent  to  subsonic  to  supersonic  components  in  the  trace 

pressure  distribution  of  a  monopole  incident  on  a  plane  as  a  function  of 
normalized  monopole  standoff  distance.  The  increase  in  standoff  causes  a 
reduction  in  evanescent  energy .  41 

3.17  Setup  and  geometry  for  the  monopole  pressure  radiation  using  the  holog¬ 
raphy  array .  43 

3.18  Setup  for  panel  excitation  from  a  monopole.  The  monopole  produces  a 

pressure  distribution  on  the  panel,  and  the  resulting  vibration  is  measured 
using  the  laser  vibrometer .  44 

3.19  Experimental  setup  for  monopole-excited  panel.  The  monopole  source  is 

the  forefront  object,  with  the  driver  (light  blue)  closest  to  the  camera.  The 
panel  is  mounted  in  the  wooden  frame;  the  white  surrounding  material  is 
Thinsulate,  which  is  wrapped  around  the  frame  to  reduce  reflections  from 
the  frame.  The  experiment  took  place  in  the  hemianechoic  chamber,  which 
reduces  reflections  from  other  sources .  45 


IX 


Figure  Page 

3.20  Wavenumber-frequency  plots  of  panel  vibration  at  kx  =  0  (left)  and  ky  =  0 

(right) .  46 

3.21  Two-dimensional  wavenumber  spectrum  of  panel  vibration  at  1000  Hz. 

Note  that  the  scales  for  the  kx  and  ky  dimensions  are  not  equal.  Note  that 

the  x  and  y  axes  are  not  equal .  46 

3.22  Monopole  pressure  amplitude  at  the  distances  measured  by  the  holography 

array.  The  1/R  decay  in  pressure  is  clearly  evident .  47 

3.23  Wavenumber  spectrum  of  the  monopole  sound  radiation  at  63  millimeter 
offset.  The  assymetric  nature  of  the  wavenumber  spectrum  is  due  to  the 
offset  measurement  of  sound  pressure  radiation.  The  black  circle  indicates 
the  material  wavenumber  in  air  of  18.3  rad/m;  the  wavenumber  spectrum 

is  clearly  highest  in  this  wavenumber  bin .  48 


3.24  Measured  ratio  of  evanescent  to  non-evanescent  components  for  the  monopole 


source  at  1000  Hertz .  49 

3.25  Wavenumber  spectrum  of  monopole  excitation  for  standoff  distances  of  22 
millimeters  (left)  and  202  millimeters  (right).  The  decreasing  energy  in 
the  system,  particularly  at  subsonic  wavenumbers  (which  is  demarcated 
by  the  dashed  black  line),  shows  the  effectiveness  of  .  49 

4.1  Schematic  of  multi-source  array  used  to  fit  pressure  distributions,  showing 

variables  Lx,  /,  9r  +  j 9t ,  Z0,  Ls,  and  N .  56 

4.2  Sample  approximate  pressure  distribution  for  a  plane  wave  pressure  dis¬ 
tribution;  approximate  pressure  distributions  with  a  fit  of  0.883  and  0.999 

are  shown .  58 

4.3  Fit  coefficient  variation  for  the  pressure  (left)  and  normal  velocity  (right) 
with  angle  and  frequency  for  a  source  array  5  meters  from  the  plane,  with 

64  sources  spanning  4  meters .  61 

4.4  Border  of  R  =  0.98  fit  plotted  with  best-fit  line  of  trace  wavelength,  Xx  = 
c/(/sin(0))  =  0.135.  The  fit  of  the  trace  wavelength  curve  shows  that  the 
trace  wavelength  is  only  an  approximate  indicator  of  the  good  fit  region.  62 

4.5  Fit  coefficient  variation  for  the  pressure  (left)  and  normal  velocity  (right) 
with  angle  and  frequency  for  a  source  array  5  meters  from  the  plane,  with 

64  sources  spanning  10  meters .  63 

4.6  Fit  coefficient  variation  for  the  pressure  (left)  and  normal  velocity  (right) 
with  angle  and  frequency  for  a  source  array  30  meters  from  the  plane,  with 

64  sources  spanning  10  meters .  63 


X 


Figure  Page 

4.7  Fit  coefficient  variation  for  the  pressure  (left)  and  normal  velocity  (right) 

with  standoff  and  source  length  for  an  incident  wave  of  2  kilohertz  at  an 
incident  angle  of  20  degrees  modeled  using  64  sources .  64 

4.8  Fit  coefficient  variation  for  the  pressure  (left)  and  normal  velocity  (right) 

with  standoff  and  source  length  for  an  incident  wave  of  5  kilohertz  at  an 
incident  angle  of  20  degrees  modeled  using  64  sources .  65 

4.9  Fit  coefficient  variation  for  the  pressure  (left)  and  normal  velocity  (right) 

with  standoff  and  source  length  for  an  incident  wave  of  5  kilohertz  at  an 
incident  angle  of  40  degrees  modeled  using  64  sources .  66 

4.10  Fit  coefficient  variation  for  the  pressure  (left)  and  normal  velocity  (right) 
with  source  length  and  number  of  sources  for  an  incident  wave  of  2  kilo¬ 
hertz  at  an  incident  angle  of  20  degrees  and  a  source  standoff  of  30  meters.  67 

4.11  Fit  coefficient  variation  for  the  pressure  (left)  and  normal  velocity  (right) 
with  source  length  and  number  of  sources  for  an  incident  wave  of  8  kilo¬ 
hertz  at  an  incident  angle  of  20  degrees  and  a  source  standoff  of  10  meters.  68 

4.12  Fit  coefficient  variation  for  the  pressure  (left)  and  normal  velocity  (right) 

with  real  and  imaginary  angle  components  for  an  incident  evanescent  wave 
of  2  kilohertz  and  a  source  geometry  of  4  meters  with  a  standoff  of  20 
meters .  69 

4.13  Fit  coefficient  improved  between  the  dipole  and  monopole  arrays  for  a 
source  array  5  meters  from  the  plane,  with  64  sources  spanning  10  meters.  70 

4.14  Fit  coefficient  improved  between  the  dipole  and  monopole  arrays  for  an 

incident  wave  of  2  kilohertz  at  an  incident  angle  of  20  degrees  and  a  source 
standoff  of  30  meters .  70 

5.1  Schematic  diagram  of  a  sample  N  —  5  multi-layer  material.  A  semi-infinite 

fluid  exists  on  each  side  of  the  three-component  material .  76 

5.2  Wave  propagation  characteristics  and  geometry,  showing  the  four  waves 
in  a  solid  system:  two  forward-going  waves,  both  longitudinal  and  shear, 

and  their  companion  backward-going  waves .  77 

5.3  Comparison  of  the  longitudinal  (left)  and  shear  (right)  transmission  coef¬ 
ficients  for  the  fluid-solid  interface  obtained  using  the  model  in  Chapter 
3  (top  row)  compared  to  that  using  the  multi-layer  model  (bottom  row). 

That  the  plots  along  each  column  match  is  indicative  of  similar  results 
between  both  models .  87 


XI 


Figure  Page 

5.4  Transmission  loss  comparison  of  Wave  Potential,  Bernoulli- Euler  and  Tim¬ 
oshenko  theories,  and  a  split  Wave  Potential  model  for  2  and  3  sublayers 
comprising  the  panel,  for  a  2cm  panel  with  30  degree  acoustic  incidence.  89 

5.5  Wavenumber-frequency  schematic  showing  phase  speed,  group  speed,  cut- 

on  frequency  relationshpis,  and  standing  waves .  91 

5.6  Representation  of  incident  angles  of  propagation  in  wavenumber-frequency 
plots.  Lines  of  equal  angle  are  shown  in  both  the  frequency-angle  plot  (left) 

and  frequency- wavenumber  plot  (right) .  92 


5.7  Geometries  of  modeled  multi-layer  systems.  The  top  system  consists  of 
10  centimeters  of  a  surrogate  layer  (and  will  also  be  used  for  the  single 
fluid  case).  The  second  system  is  the  same  surrogate  layer  enclosed  by 
thin  polycarbonate  walls,  to  simulate  a  container.  The  lowest  geometry  is 
the  same  surrogate  and  walls,  only  with  a  1  centimeter  air  gap  between 


the  surrogate  and  the  polycarbonate  wall .  95 

5.8  Intensity  in  middle  layer  of  the  fluid  layer  system,  with  the  wave  speed  of 
the  fluid  (700  meters  per  second)  and  the  wave  speed  in  air  (343  meters 

per  second)  indicated .  96 

5.9  Intensity  in  the  middle  layer  of  the  solid  layer  system,  with  rclevent  fea¬ 
tures  and  wave  speeds  labeled .  97 

5.10  Intensity  in  the  center  of  the  middle  (surrogate)  layer  of  the  plastic- 

bounded  surrogate  system .  98 

5.11  Intensity  in  the  middle  layer  of  middle  (surrogate)  layer  of  the  plastic- 

bounded  system  with  an  air  gap  between  the  plastic  and  surrogate.  .  .  99 

5.12  Intensity  transmission  loss  in  the  center  layer  (left)  and  transmission  loss 

improvement  due  to  the  evanescent  wave  (right)  for  the  single-layer  sur¬ 
rogate  layer  system  with  0.1  radian  evanescent  component .  100 

5.13  Intensity  transmission  loss  for  the  polycarbonate-shielded  surrogate  sys¬ 
tem,  consisting  of  polycarbonate  in  front  of  10  centimeters  of  surrogate. 
Intensity  transmission  loss  is  shown  in  the  polycarbonate  layer  (left)  and 

the  surrogate  layer  (right) .  101 


5.14  Intensity  transmission  loss  with  evanescent  component  for  the  polycarbonate- 


shielded  surrogate  system,  consisting  of  polycarbonate  in  front  of  10  cen¬ 
timeters  of  surrogate.  Intensity  transmission  loss  is  shown  in  the  polycar¬ 
bonate  layer  (left)  and  the  surrogate  layer  (right) .  102 

5.15  Intensity  transmission  reduction  in  the  surrogate  layer  for  the  polycarbonate- 

shielded  surrogate  system .  102 


Figure  Page 

6.1  Surrogate  samples  tested.  The  50%  sample  is  on  top,  consisting  of  lower 
stiffness  and  a  more  clastic  structure.  The  75%  sample  is  significantly 
more  rigid,  and  would  fracture  under  substantial  deformation .  107 


6.2  Schematic  of  experimental  measurements  of  polycarbonate  beam.  Two 

excitation  locations  were  used,  with  the  ’’side”  excitation  being  along  the 
same  plane  as  the  measurement,  and  the  ’’end”  excitation  being  at  the 
far  end  of  the  polycarbonate  bar.  For  mounting,  the  beam  was  hung 
(suspended  from  the  end  above  a  fixed  point),  given  a  foam  backing  that 
covered  the  entirety  of  the  surface  opposite  the  measurement,  and  posi¬ 
tioned  with  minimal  support  that  allowed  for  reduced  contact  with  the 
sample .  109 

6.3  Polycarbonate  vibration  for  the  hung  (top  row),  foam  backing  (middle 
row)  and  minimal  foam  support  (bottom  row)  conditions,  with  side  (left 
plots)  and  end  (right  plots)  shown.  Significant  features  have  been  labeled 


on  the  plot .  113 

6.4  50  %  surrogate  vibration  for  foam  backing  (top)  and  polycarbonate  back¬ 

ing  (bottom)  conditions,  with  side  (left  plots)  and  end  (right  plots)  shown. 
Significant  features  have  been  labeled  on  the  plot .  115 

6.5  Vibration  in  the  position-frequency  domain  for  the  50%  surrogate  sample, 

comparing  the  foam  mounting  condition  for  side  (left)  and  end  (right) 
excitation .  116 

6.6  75  %  surrogate  vibration  for  foam  backing  (top)  and  polycarbonate  back¬ 

ing  (bottom)  conditions,  with  side  (left  plots)  and  end  (right  plots)  shown. 
Significant  features  have  been  labeled  on  the  plot .  117 


Xlll 


ABSTRACT 

Jessop,  Andrew  M.  Ph.D,  Purdue  University,  December  2013.  Near-Field  Pressure 
Distributions  to  Enhance  Sound  Transmission  into  Multi-Layer  Materials.  Major 
Professor:  J.  Stuart  Bolton,  School  of  Mechanical  Engineering. 

The  large  impedance  difference  between  air  and  most  solids  prevents  significant 
energy  transfer  from  incident  acoustic  waves  across  the  air-material  interface.  Refrac¬ 
tion  also  plays  a  role  in  preventing  acoustic  transmission,  as  the  wave  speed  difference 
between  air  and  solid  materials  results  in  an  increase  of  the  resulting  propagation 
angles,  creating  near-field  pressure  distributions  in  the  solid  material.  By  utilizing 
evanescent  pressure  distributions,  which  decay  normal  to  the  usual  direction  of  prop¬ 
agation  and  are  represented  as  plane  waves  propagating  with  complex  angles,  energy 
propagation  through  the  interface  can  be  increased  in  the  subsonic  region  of  wave 
propagation:  i.e.,  where  waves  typically  do  not  propagate  into  a  material  with  any 
effectiveness.  By  using  an  array  of  sources,  it  is  possible  to  produce  evanescent  pres¬ 
sure  distributions  in  the  solid.  The  way  in  which  the  characteristics  of  this  array  of 
sources  affect  the  efficiency  of  the  generation  of  evanescent  pressure  distributions  are 
explored. 

Because  high  impedance  materials  can  be  paired  to  a  lower  impedance  materials 
of  interest  to  impede  acoustical  energy  transmission,  the  wave  propagation  through 
multi-layer  materials  must  be  considered  to  give  an  full  accounting  of  power  transmis¬ 
sion  into  structures.  A  model  for  wave  propagation  in  multi-layer  systems  of  solids  and 
fluids  was  developed  using  wave  potentials  in  each  layer,  allowing  for  coupling  between 
material  types  and  calculation  of  inter-layer  states.  By  using  wavenumber-frequency 
analysis,  it  is  possible  to  target  specific  components  in  the  multi-layer  system  and 
understand  the  particular  wavetypes  that  cause  energy  propagation  into  the  system. 
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CHAPTER  1.  INTRODUCTION 

Evanescent  waves,  i.e.,  waves  that  simultaneously  propagate  and  attenuate,  can  be 
used  to  improve  energy  transmission  into  high  impedance  solid  materials.  By  over¬ 
coming  restrictions  due  to  refraction  at  the  air-solid  interface,  an  evanescent  wave 
incident  on  the  surface  can  propagate  energy  into  a  solid  material  at  a  greater  range 
of  angles  than  classical  plane  waves.  For  high-impedance  materials  in  multi-layer 
systems,  geometry  effects  may  provide  addditional  avenues  for  energy  transmission 
into  the  material.  An  understanding  of  evanescent  wave  effects  and  their  propagation 
into  multi-layer  materials  will  enable  novel  methods  for  excitation  of  those  materials. 

Acoustical  transmission  into  solid  materials  is  typically  limited  by  the  large  impedance 
difference  between  air  and  the  material.  Lower-density  solid  materials  have  densities 
at  least  1,000  times  greater  than  air,  and  speed  of  sound  differences  of  10  times  or 
more  are  typical.  While  the  density  difference  can  be  overcome  by  exploiting  the 
compressional  and  shear  coupling  effects  in  the  material,  refraction  due  to  the  differ¬ 
ence  in  sound  speed  is  a  limit  that  prevents  energy  transmission  beyond  the  critical 
angle  of  incidence. 

The  production  of  evanescent  pressure  distributions  requires  special  sources.  When 
using  a  multi-source  array,  each  source  can  be  given  a  phase  and  amplitude  in  such  a 
way  that  the  interference  of  the  wavefronts  produces  the  desired  wave  characteristics 
on  a  surface.  The  amplitudes  and  phases  of  these  sources  can  be  found  by  construct¬ 
ing  a  least-squares  calculation  of  source  strength  that  would  approximate  a  wave  at  a 
particular  angle  of  incidence.  Factors  such  as  the  angle  and  frequency  of  the  desired 
pressure  distribution,  as  well  as  the  geometry  of  the  sources,  all  affect  the  ability  of 
the  source  array  to  produce  an  approximation  of  the  desired  pressure  distribution. 

For  more  complex  multi-layer  structures,  the  geometry  and  component  materials 
both  affect  the  wave  transmission  into  the  system.  Accounting  for  geometry  will  allow 


2 


for  flexural  waves  and  cut-on  of  modes  in  each  layer,  which  will  affect  the  total  sound 
transmission  through  the  system. 

The  wave  propagation  characteristics  of  materials  with  unknown  properties  can 
be  found  experimentally  and  used  to  improve  the  accuracy  of  wave  propagation  mod¬ 
els.  The  waves  resulting  from  such  tests  vary  greatly  with  excitation  location,  and 
the  materials  tested  need  to  be  mounted  properly  to  prevent  damping  and  geom¬ 
etry  constraints.  Understanding  these  effects  and  the  waves  produced  by  different 
experimental  conditions  will  allow  for  assessment  of  wave  propagation  and  material 
properties. 

1.1  Motivation 

The  research  undertaken  in  this  project  was  part  of  a  multi-university  research 
initiative  (MURI)  focused  on  non-contact  excitation  of  energetic  materials.  There  are 
existing  spectrographic  methods  for  large-standoff  detection  of  chemical  signatures, 
but  when  these  chemical  signatures  are  weak,  excitation  is  needed  to  produce  more 
outgassing  of  the  energetic  material.  Electromagnetic  methods  of  excitation  can  be 
defeated  be  installing  a  dielectric  buffer  around  the  material;  in  those  cases,  to  excite 
the  materials,  mechanical  displacement  is  necessary.  Acoustical  excitation  may  be 
able  to  provide  a  large-standoff  method  for  excitation. 

A  typical  scenario  involves  an  energetic  material  behind  a  thin,  rigid  material. 
Such  a  material  would  serve  as  both  a  container  for  the  energetic  material  (which 
typically  needs  to  be  cast  from  a  liquid  state)  as  well  as  concealment  of  the  energetic 
material.  Examples  would  be  a  milk  jug,  which  is  composed  of  a  thin  layer  of  plastic, 
or  an  aluminum  can.  There  may  be  additional  barriers  serving  as  concealment,  as 
well  as  imperfections  in  the  coupling  between  the  barrier  materials  and  the  energetic 
material  that  produce  air  gaps  and  further  limit  energy  transmission. 

Related  aspects  of  the  research  performed  by  other  reseachers  in  this  overall 
project  involve  classification  of  the  coupling  of  the  materials,  its  response  to  low 
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frequency  excitations,  and  extensive  exploration  of  the  thermo-physical  properties 
that  may  allow  for  heating  and  outgassing. 

1.2  Contribution  of  Research 

The  primary  goal  of  this  research  was  the  exploration  of  evanescent  wave  energy 
transfer  into  multi-layer  materials.  To  that  end,  several  methodologies  to  increase  our 
understanding  of  evanescent  waves  and  techniques  for  simulation  of  wave  propagation 
through  multi-layer  materials  have  been  developed. 

A  detailed  discussion  of  evanescent  waves  in  air,  as  a  result  of  complex  angle 
of  transmission  along  a  plane,  will  allow  for  the  development  of  a  model  exploring 
energy  transfer  of  evanescent  waves  across  the  fluid-solid  interface.  While  this  model 
is  based  on  previous  formulations  of  fluid-solid  interaction,  an  accounting  for  the 
effects  of  evanescent  waves  has  not  previously  been  undertaken. 

Wave  propagation  in  multi-layer  materials  is  modeled  by  using  wave  potentials  in 
each  layer  to  allow  for  inter-layer  material  state  calculation.  In  the  model,  the  incom¬ 
ing  and  outgoing  waves  are  separated,  which  allows  for  calculation  of  the  intensity  of 
each  wave  and  a  more  accurate  understanding  of  the  effects  of  each  layer  on  energy 
propagation.  In  addition,  since  the  incoming  and  outgoing  waves  are  separated,  an 
evanescent  incident  wave  can  be  accurately  represented. 

A  detailed  accounting  is  also  given  of  the  effects  of  excitation  and  material  mount¬ 
ing  conditions  on  wave  propagation  measurements.  This  information  can  be  used  to 
improve  classification  of  wavetypes  in  unknown  materials. 

Generation  of  plane  and  evanescent  waves  using  multi-source  methods  was  also 
explored  and  a  simple  model  was  developed.  The  contribution  of  the  research  on  this 
subject  is  a  thorough  exploration  of  the  effects  of  the  incident  wave  properties  and 
source  geometry  on  the  ability  to  produce  desired  pressure  distributions  in  energetic 
materials. 
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1.3  Outline 

After  a  discussion  of  previous  work  relevent  to  this  research  is  given  in  Chapter 
2,  Chapter  3  begins  with  a  discussion  of  evanescent  waves  and  how  they  can  be 
represented  using  complex  angles.  An  evanescent  wave’s  ability  to  increase  energy 
transmission  in  the  subsonic  region  of  incident  angles  will  be  presented  using  closed- 
form  solutions  of  fluid-solid  interaction.  An  assessment  of  the  evanescent  waves  that 
are  generated  in  the  near-field  of  a  monopole  then  allows  for  experimental  exploration 
of  their  effects  on  energy  transmission  into  solids. 

Chapter  4  is  focused  on  multi-source  methods  for  generating  evanescent  pressure 
distributions.  It  will  begin  with  a  discussion  of  the  radiation  effects  of  simple  sources 
and  the  incident  wave  that  is  to  be  modeled,  before  presenting  a  least-squares  formu¬ 
lation  for  approximating  incident  waves  using  multiple  simple  sources.  An  analysis 
of  the  effects  of  incident  wave  parameters  and  source  geometry  follows,  including  a 
comparison  of  monopole  and  dipole  efficiency  in  modeling. 

In  Chapter  5  a  model  for  multi-layer  wave  propagation  using  wave  potentials 
is  presented.  Several  representative  models  are  explored  to  show  how  multi-layer 
systems  affect  wave  propagation. 

Chapter  6  contains  a  description  of  experimental  tests  of  wave  propagation  in  both 
previously  classified  and  unknown  materials.  These  tests  allowed  for  an  exploration  of 
mounting  and  excitation  effects  on  the  resulting  vibration  in  the  material,  and  their 
effects  on  wave  properties. 

Finally,  a  summary  of  the  research  findings  and  recommendations  for  continuing 
avenues  of  research  are  given  in  Chapter  7. 
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CHAPTER  2.  PREVIOUS  RESEARCH 

The  results  of  the  literature  survey  are  organized  into  three  general  areas  of  research: 
near-field  effects  on  sound  transmission;  wave  propagation  in  multi-layer  systems;  and 
generation  of  arbitrary  pressure  distributions.  The  chapter  ends  with  a  summary  of 
what  is  and  is  not  known  is  these  areas  that  is  relcvent  to  the  current  research. 

2.1  Near-Field  Effects  on  Sound  Transmission 

Much  research  has  been  performed  regarding  acoustic  transmission  into  solids  us¬ 
ing  waves  in  the  ultrasonic  frequency  domain.  For  example,  Gan,  Hutchins,  Billson, 
and  Schindel  [1]  have  looked  at  the  transmissibility  of  ultrasonic  waves  through  ma¬ 
terial  samples.  Many  of  these  methods  exploit  resonances  of  the  structure  or  the 
microstructure  of  the  material  being  investigated.  However,  ultrasound  is  only  useful 
in  applications  where  propagation  distances  through  air  are  short  because  the  decay 
of  ultrasonic  waves  in  air  is  significant;  Fox,  Khuri- Yakub,  an  Kino  [2],  for  example, 
found  decay  levels  of  52  dB  over  40  centimeters.  Because  a  primary  goal  of  the  work 
reported  here  is  large-standoff  investigation  of  materials,  ultrasonic  methods  will  not 
be  useful. 

The  basic  theory  of  solid-fluid  interaction  is  explained  by  Brekhovskikh  [3].  By 
equating  the  pressure  and  velocity  in  the  wave-carrying  fluid  with  that  on  the  solid 
surface,  Brekhoskikh  derived  the  equations  for  acoustic  reflection  and  transmission 
that  are  used  as  the  basis  for  this  study.  He  derived  transmission  coefficients  for  both 
the  compressional  and  shear  waves  propagating  into  a  solid,  and  performed  a  thorough 
exploration  of  various  regions  of  transmission,  such  as  the  supersonic  and  subsonic 
regions  of  both  the  compressional  and  shear  waves.  Brekhovskikh  also  derived  transfer 
matrices  for  multi-layer  propagation. 
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The  most  common  high- impedance-difference  interface  studied  is  the  air-water 
interface,  due  to  its  naval  applications.  Hudimac  [4]  looked  at  the  effects  of  refraction 
across  the  interface  on  energy  transmission,  deriving  iso-intensity  lines  for  sound 
originating  in  air  and  transmitted  into  water.  Weinstein  and  Henney  [5]  used  that 
work  as  the  basis  to  calculate  monopole  sound  transmission  into  water,  looking  at 
energy  transmission  as  a  function  of  height  above  the  surface.  Young  performed 
experimental  measurements  of  sound  transmission  between  water  and  air  [6].  From 
his  measurements,  he  found  pressure  levels  equal  to  those  predicted  by  a  theory  based 
on  ray-tracing,  with  an  accuracy  of  1-2  dB  at  most  levels. 

Meecham  looked  at  sound  propagation  into  water  from  an  airborne  source  [7]. 
He  formulated  conditions  in  which  acoustical  energy  from  the  air  would  propagate 
into  the  water  at  small  grazing  angles  caused  by  refraction  in  the  air.  He  also  found 
transmission  enhancement  effects  due  to  surface  roughness. 

Bobrovnitskii  [8,  9]  briefly  discussed  the  possibility  of  energy  flow  from  evanescent 
waves.  He  shows  that  the  power  resulting  from  a  combination  of  evanescent  waves 
in  opposing  directions  will  lead  to  interaction  between  the  otherwise  out-of-phase 
velocity  and  pressure  components,  resulting  in  an  intensity  equal  to  the  cross-product 
of  the  two  waves’  amplitudes. 

Godin  [10,  11,  12,  13]  has  explored  theoretical  conditions  under  which  a  high- 
impedance  difference  boundary  becomes  acoustically  transparent  due  to  evanescent 
waves  from  a  source.  In  his  research  on  both  the  air-water  interface  [10,  11,  12]  and 
the  air-solid  interface  [13],  he  formulated  equations  for  transmission  from  a  monopole 
source  in  the  higher-impedance  solid  material  into  the  fluid.  Godin  found  significant 
power  flux  across  the  interface  due  to  inhomogeneous  wave  components  in  a  monopole 
close  to  the  interface.  For  the  air-solid  transmission  interface,  Godin  also  accounted 
for  various  lossy  wavetypes  that  are  typically  seen  in  this  type  of  transmission.  He  also 
investigated  the  effects  of  interface  roughness;  he  found  that  there  is  no  significant 
additional  transmission  due  to  roughness  that  is  much  smaller  than  a  wavelength  in 
the  incident  fluid  [10].  Godin’s  theory  was  tested  by  Calvo,  McDonald,  Nicholas, 
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and  Orris,  who  further  developed  the  formulation  to  account  for  total  power  flux 
when  accounting  for  reflections  from  the  surface  [14],  and  performed  preliminary 
experimental  measurements  [15].  However,  in  these  formulations  the  authors  made 
use  of  inhomogeneous  wave  components  in  water  which  are  supersonic  in  air;  the 
inverse  of  this  problem  cannot  take  advantage  of  such  a  range  of  inhomogeneous 
waves. 

The  coupling  of  acoustical  energy  and  flexural  waves  in  materials  has  been  de¬ 
scribed  by  Fairy  and  Gardonio  [16].  Flexural  waves  typically  have  much  lower  wave 
speeds  than  longitudinal  or  shear  waves  in  a  material,  and  they  can  therefore  be 
excited  across  a  much  broader  range  of  incident  angles.  Flexural  waves  are  not  use¬ 
ful  for  formulation  of  sound  power  transmission  across  a  material  interface,  however, 
because  the  flexural  wave  is  structure-dependent  and  will  not  provide  general  trans¬ 
mission  across  a  variety  of  different  geometries. 

Park  et  al.  [17]  discussed  the  use  of  metamaterials  in  increasing  the  transfer  of 
electromagnetic  waves  into  materials.  By  using  a  structure  with  an  effective  cell  size 
much  smaller  than  a  wavelength,  a  negative  density  that  enhances  the  transmission 
of  an  evanescent  electromagnetic  wave  can  be  created.  However,  despite  the  use  of 
metamaterials  being  shown  to  increase  the  acoustic  transmission  in  classical  wave 
propagation  [18],  these  evanescent  transmission  effects  cannot  be  seen  in  acoustical 
materials  because  mechanical  dissapation  negates  the  effects  of  the  negative  density 
in  the  metamaterial. 

2.2  Wave  Propagation  in  Multi-Layer  Systems 

Folds  and  Loggins  [19]  expanded  on  the  sound  transmission  work  performed  by 
Brekhovskikh  [3]  to  create  propagation  matrices  for  general  cases  of  multi-layer  solids. 
He  calculated  ultrasonic  wave  propagation  through  single-layer  and  three-layer  panels, 
and  found  good  correlation  between  the  model  and  experimental  results  for  lower 
grazing  angles. 


Pierce  [20]  explains  more  general  theories  for  multi-layer  propagation,  such  as  the 
effects  of  inter-layer  states  and  the  ideal  properties  for  a  coupling  layer.  He  also 
considers  flexural  waves,  though  as  an  independent  wavetype  rather  than  through 
coupling  effects  of  material  waves  in  a  thin  plate  system. 

Brouard,  Lafarge,  and  Allard  formulated  a  system  of  linear  equations  to  model 
sound  propagation  in  mixed  fluid-solid  layered  systems  [21].  By  working  from  similar 
matrix-based  methods  for  sound  propagation  in  fluids  and  solids,  Brouard  developed 
a  matrix  relating  the  pressure  and  velocity  of  waves  propagating  through  multiple 
layers.  At  interfaces  where  the  medium  changes,  such  as  from  fluid  to  solid,  a  coupling 
matrix  is  used  to  equate  the  different  types  of  pressure  and  velocity  amplitudes. 
The  use  of  this  formulation  allows  for  generation  of  transfer  functions  between  the 
incoming  pressure  and  the  resulting  pressure  or  velocity  at  an  internal  layer. 

Auld’s  investigated  propagation  in  waveguides,  and  derived  dispersion  relation¬ 
ships  for  the  different  wavetypes  in  the  system  [22],  In  his  calculations  using  exist¬ 
ing  models,  based  primarily  on  Folds  and  Pierce’s  formulations,  he  found  dispersion 
relationships  showing  material  and  flexural  waves  in  constrained  systems.  He  ex¬ 
plained  the  relationship  of  cut-on  frequencies  for  particular  modes  of  vibration  and 
how  they  transitioned  into  waveguide-like  behavior  at  higher  angles  of  incidence.  He 
also  showed  how  coupling  of  material  waves  can  create  Lamb  waves,  such  as  flexural 
and  dilatational  waves. 

Maidanik  and  Dickey  developed  a  methodology  for  deriving  a  low-reflectivity  ma¬ 
terial  [23].  In  a  two-layer  material  consisting  of  a  thin  plate  and  a  thicker  mate¬ 
rial  layer  underneath,  certain  conditions  for  resonance  could  be  met  that  couple  the 
flexural  waves  in  the  plate  with  other  wavetypes  in  the  material  to  produce  a  low- 
reflectivity  condition.  They  explored  the  material  combinations  under  which  this 
could  occur  as  well. 

Mead  has  looked  at  vibration  propagation  through  periodic  systems  [24,  25].  He 
accounted  for  the  presence  of  the  modes  of  individual  layers  of  a  multi-layer  sys- 
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tem,  and  found  resonance  conditions  that  would  cause  ideal  wave  propagation  in  the 
system. 

In  previous  work  by  the  author  [26]  on  assessing  the  effects  of  deformation  on  the 
tire  acoustic  mode,  the  author  used  a  wave  potential-based  method  in  a  system  of 
linear  equations.  The  work  was  built  upon  a  model  developed  by  Thompson  [27],  who 
used  wave  potentials  across  a  simpler  geometry  to  formulate  a  closed-form  solution 
for  the  effects  of  tire  deformation  on  the  tire’s  acoustic  mode. 

Grosh  and  Williams  developed  a  wavenumber-frequency  decomposition  technique 
to  visualize  wavetypes  in  a  material  [28].  In  addition  to  other  techniques  such  as 
Prony  series  decomposition,  often  used  in  modal  analysis  of  vibrations  in  structures, 
the  use  of  wavenumber-frequency  decomposition  allows  for  identification  of  vibrational 
modes.  Grosh  et  al.  applied  this  technique  to  both  measured  and  simulated  data. 

Wahl  and  Bolton  calculated  wave  speed  properties  for  dispersive  and  non-dispersive 
waves  in  a  system  using  Wigner  distributions  [29].  By  visualizing  the  amplitude  of  a 
pulse  versus  time  and  frequency,  the  arrival  time  of  waves  and  the  resulting  speed  of 
wave  propagation  can  be  determined.  By  filtering  elements  of  the  Wigner  distribu¬ 
tion,  they  were  able  to  determine  the  energy  carried  by  various  wavetypes. 

Bolton,  Kim,  and  Song  others  have  used  a  wavenumber  domain  representation  of 
vibration  to  visualize  wavetypes  in  a  system  [30,  31,  32],  By  visualizing  tire  vibra¬ 
tion  in  the  wavenumber-frequency  domain,  the  wavetypes  specific  to  tire  treadband 
vibration  were  isolated,  and  the  damping  in  the  treadband  could  be  assessed. 

In  the  previous  research  discussed  here,  phenomena  evident  in  anomalous  fluid- 
solid  transmission  has  been  identified.  By  further  classifying  these  effects,  there  ex¬ 
ists  the  possibility  of  controlling  and  utilizing  them  to  deliver  desired  excitations  to 
a  multi-layer  material.  The  combination  of  the  previously-developed  methodologies 
for  simulating  multi-layer  acoustical  propagation  with  the  technique  of  wavenumber- 
frequency  visualization  will  allow  for  isolation  of  particular  wavetypes  in  the  system. 
Incident  evanescent  can  be  created  by  expanding  upon  the  techniques  shown  in  pre¬ 
vious  research  of  wave  front  synthesis,  considered  next. 
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2.3  Generation  of  Arbitrary  Pressure  Distributions 

Typical  multi-source  methods  entail  the  use  of  an  array  of  sources  with  selectively- 
phased  outputs  designed  to  produce  an  optimal  amplitude  output.  For  example, 
beamforming  techniques  [33,  34]  are  applied  to  source  arrays  with  the  goal  of  ampli¬ 
fying  an  output  signal  in  a  direction  of  interest.  However,  the  present  focus  was  on 
creating  a  specified  pressure  distribution,  rather  than  amplifying  a  source. 

In  other  multi-source  methods  the  array  of  sources  was  often  used  to  generate 
a  plane  wave:  i.e.,  the  sources  were  phased  to  generate  a  relatively  large  region  of 
plane  wavefront.  Berkhout,  De  Vries,  and  Vogel  looked  at  using  multiple  sources  to 
generate  plane  waves  in  the  far-held  [35]  by  using  a  loudspeaker  array.  By  selectively 
phasing  the  sources,  a  larger  wavefront  could  be  generated  by  combining  the  smaller 
wavefronts  of  the  component  sources.  Their  research  was  primarily  focused  on  sim¬ 
ulating  large  wavefronts  at  a  particular  location;  such  a  method  could  not  generally 
produce  plane  waves  at  an  angle  to  an  interface  surface. 

Chang,  Choi  and  Kim  looked  at  generating  plane  waves  in  free  space  using  a  spher¬ 
ical  array  of  loudspeakers  [36].  In  his  method,  he  generated  a  plane  wave  held  using 
multiple  sound  sources,  but  the  dependence  on  the  spherical  geometry  surrounding 
the  region  of  interest  makes  it  unfeasible  for  our  purposes. 

Russell,  Titlow  and  Bemmen  [37]  looked  at  the  power  output  and  directivity 
patterns  for  monopoles,  dipoles,  and  quadrupoles  and  measured  how  the  theoretical 
formulations  of  such  sources  match  experimental  approximations.  He  found  that  his 
experimental  approximations  matched  theory  reasonably  well  at  low  frequencies,  but 
were  limited  by  geometry  and  diffraction  at  higher  frequencies. 

Bolton,  Beauvilain,  and  Gardner  experimented  with  the  generation  of  pressure  dis¬ 
tributions  in  a  region  by  using  simple  sources,  particularly  dipoles  and  quadrupoles 
[38].  Through  selective  phasing  of  colocated  simple  sources,  the  radiation  charac¬ 
teristics  of  the  system  could  be  varied,  including  variation  of  the  angular  pressure 
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distribution  and  reduction  of  the  near-field  intensity.  These  theories  were  verified  in 
later  work  [39]  using  four  uniquely-arranged  sources. 

Ahrens  and  Spors  developed  methods  for  reproduction  of  plane  waves  using  circu¬ 
lar  [40]  and  linear  arrays  [41].  In  their  formulation,  they  discretized  continuous  sound 
sources  to  generate  the  desired  plane  wave  field.  By  calculating  a  driving  function  for 
the  continuous  distribution  that  gave  an  exact  answer  for  the  pressure  distribution, 
then  discretizing  that  driving  function  into  an  array  of  point  sources,  they  were  able 
to  approximate  the  sound  field  produced  by  the  discrete  source.  They  also  studied 
driving  function  and  discretization  techniques  to  minimize  the  error  associated  with 
their  linearization  method.  These  theoretical  formulations  will  be  difficult  to  adapt 
to  loud  speaker  arrays  used  in  experimental  testing  due  to  the  complexity  of  the  driv¬ 
ing  functions;  a  methodology  that  allows  for  arbitrary  angular  variation  of  a  discrete 
source  would  be  more  useful. 

Robin,  Berry,  Moreau,  and  Dia  [42]  used  multi-source  methods  to  generate  repro¬ 
ducible  turbulent  pressure  fields  for  use  in  structural  acoustics.  In  their  approach  they 
used  both  wavefront  synthesis  (as  defined  by  Berkhout)  and  near-field  holography  the¬ 
ory  (inverted  from  its  typical  use  in  sound  source  detection)  to  generate  the  wave  field 
on  a  surface  by  using  controllable  sources.  To  simulate  the  turbulent  boundary  layer 
properties  by  using  a  limited  array  of  variables,  they  used  relationships  between  the 
desired  wavelength  of  reproduction  and  the  number  of  sources  used. 

Wang  and  Wu  formulated  a  method  for  generating  pressure  distributions  resulting 
from  irregularly-shaped  vibrating  sources  [43].  They  generated  a  set  of  orthonormal 
functions  from  the  modes  of  structure,  then  formed  those  equations  into  a  least- 
squares  solution  to  calculate  the  combination  of  modes  that  produced  the  measured 
pressure  distribution.  Wu  expanded  on  this  approach  and  applied  it  to  a  radiating 
structure  in  an  enclosed  vibrating  space,  and  combined  those  principles  with  nearfield 
acoustical  holography  to  create  a  computationally  efficient  method  for  sound  field 
prediction  [44,  45]. 
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Kirkeby  and  Nelson  formulated  a  method  of  generating  a  complex  pressure  dis¬ 
tribution  from  a  series  of  simple  sources  [46].  Their  method  consists  of  altering  the 
amplitude  and  phase  of  a  system  of  monopoles  in  a  given  geometry.  They  accom¬ 
plished  this  by  generating  a  series  of  linear  equations,  consisting  of  a  matrix  of  Green’s 
functions  between  the  sources  and  pressure  distribution  at  a  point,  multiplied  by  the 
unknown  phase  amplitudes  of  the  sources,  to  equal  the  resulting  pressure  distribution 
at  a  series  of  discrete  points.  By  finding  a  least-squares  solution  to  this  matrix,  they 
were  able  to  find  a  vector  of  source  phases  that  approximated  the  pressure  distribu¬ 
tion.  They  also  established  error  bounds  for  the  solution  for  different  geometry  and 
frequency  types. 

2.4  Conclusions 

The  work  described  in  this  thesis  on  evanescent  wave  energy  propagation  into 
solid  materials  is  an  expansion  of  the  theory  developed  by  Brekhovskikh  [3] ;  his  work 
allows  for  evanescent  incident  waves  as  part  of  the  general  solution  for  fluid-solid 
transmission  but  does  not  consider  the  consequences  of  using  such  waves.  By  formu¬ 
lating  evanescent  waves  incident  on  the  surface  by  using  complex  angles,  it  can  be 
shown  that  significant  energy  transmission  across  the  fluid-solid  interface  is  possible. 
Similiar  propagation  characteristics  were  shown  in  the  near-field  of  a  monopole  by 
Godin  [10,  11,  12,  13];  however,  in  the  present  work  further  classification  is  provided 
by  using  specifically  generated  evanescent  waves  rather  than  those  already  present  in 
the  near-field.  In  his  research  Godin  also  established  the  basis  for  the  experimental 
portion  of  the  research  (reported  here)  on  evanescent  wave  transmission  into  solid 
matrials. 

The  formulation  of  evanescent  wave  pressure  distributions  described  in  Chapter 
4  is  an  extension  of  the  work  of  Kirkeby  et  al  [46].  By  using  their  system  of  lin¬ 
ear  equations,  an  accounting  for  the  effects  of  source  geometry  and  desired  pressure 
distribution  can  be  made.  Other  work,  such  as  that  performed  by  Russell  [37]  on 
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experimental  verification  of  monopole  and  dipole  properties,  were  useful  in  assuring 
experimental  reproductivity  of  the  distribution  produced. 

The  model  for  wave  propagation  through  layers  used  in  this  research  was  initially 
developed  by  Brouard  [21],  who  expanded  the  transfer  matrix  approach  for  wave 
propagation.  His  model,  a  system  of  linear  equations,  accounts  for  both  forward-  and 
backward-propagating  waves  in  the  system.  The  use  of  wave  potentials  as  inputs  to 
the  system  of  equations  is  based  on  previous  work  by  Jessop  [26]  (which  itself  is  taken 
from  the  work  of  Thompson  [27])  to  account  for  wave  potentials  rather  than  material 
states,  thus  allowing  for  inter-layer  state  and  intensity  calculation. 
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CHAPTER  3.  NEAR-FIELD  EFFECTS  ON  SOLID-FLUID 
ACOUSTICAL  TRANSPARENCY 

The  characteristics  of  evanescent  pressure  distributions  can  be  used  to  enhance  sound 
propagation  through  the  air-solid  interface.  In  typical  air-solid  propagation,  the  dif¬ 
ference  in  density  and  wave  speed  across  the  interface  results  in  very  low  energy 
transmission.  The  difference  in  wave  speed  also  causes  extreme  refraction  through 
the  interface;  pressure  in  the  solid  will  then  decay  quickly  and  there  will  be  little 
velocity  propagation  into  the  material.  The  region  wherein  refraction  reduces  energy 
propagation  into  the  material  is  known  as  the  subsonic  region;  it  is  characterized  by 
the  normal  speed  of  sound  of  a  wave  into  the  material  being  less  than  the  material’s 
wave  propagation  speed.  By  inverting  the  typical  near-field  effects  seen  in  the  sub¬ 
sonic  region  of  acoustical  transmission  and  specifically  generating  evanescent  pressure 
distributions,  it  will  be  a  shown  that  a  larger  portion  of  the  energy  in  the  sound  wave 
can  be  transmitted  into  the  solid  material. 

Evanescent  pressure  distributions  can  be  modeled  by  a  plane  wave  propagating 
at  a  complex  angle.  As  shown  by  Brekhovskikh  [3],  the  complex  angle  represents 
decay  or  growth  in  the  wave  amplitude  normal  to  the  direction  of  propagation.  Such 
distributions  are  typically  seen  as  the  result  of  subsonic  acoustic  transmission  (i.e., 
at  incidence  angles  larger  than  the  critical  angle);  however,  intentional  generation  of 
such  a  distribution  may  provide  a  means  of  controlling  the  angle  of  propagation  and 
energy  transfer  across  an  interface. 

This  chapter  begins  with  a  description  of  evanescent  pressure  distributions  made 
by  the  inclusion  of  complex  wavenumbers  and  angles  into  the  standard  plane  wave 
formulas.  A  discussion  of  the  theory  of  acoustical  transmission  through  an  interface 
will  lead  to  an  exploration  of  the  effects  of  evanescent  pressure  distributions  in  acous- 
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tical  transmission.  The  evanescent  pressure  distributions  generated  by  a  monopole 
will  also  be  discussed  and  verified  experimentally. 


3.1  Evanescent  Plane  Waves 

3.1.1  Fundamentals  of  Wave  Propagation 

The  pressure  of  a  propagating  plane  wave  at  a  point  in  space,  x,  is: 


p  =  Aejlse~jut  ,  (3.1) 

where  A  is  the  amplitude  (and  initial  phase)  of  the  wave,  k  is  the  wavenumber  vector 
(the  rate  of  change  in  phase  per  unit  distance),  and  u  is  the  radian  frequency  (the 
rate  of  change  in  phase  over  time,  equal  to  2n  times  the  frequency  /  in  radians 
per  second).  The  components  of  k  specify  a  direction  of  propagation  for  the  wave. 
For  steady-state  formulations,  the  phase  change  with  time,  e does  not  have  an 
impact  on  the  parameters  discussed  here,  and  is  excluded  from  all  future  equations 
for  brevity. 

The  material’s  wavenumber  k  =  \k\,  which  is  the  phase  change  per  unit  distance  in 
the  direction  of  wave  propagation,  is  a  unit  inversely  proportional  to  the  wavelength 
of  sound,  A,  and  is  related  to  frequency  by  the  speed  of  sound  in  the  material  c: 


k  2tt  uj 

X  c 


(3.2) 


A  wave  propagating  over  a  plane  surface  will  have  a  trace  wavenumber  (and  wave¬ 
length)  along  the  surface  that  is  different  from  the  material  wavenumber.  If  a  wave  is 
incident  on  a  surface  oriented  along  the  x-y  plane  at  an  angle  6  from  the  normal,  z,  as 
depicted  in  Figure  3.1,  the  effective  wavenumber  along  the  surface  in  the  i'-direction 


is: 


kx  =  k  sin  6  . 


(3.3) 
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At  normal  incidence,  kx  is  0,  representing  a  uniform  phase  over  the  entire  incident 
plane.  With  increasing  angle,  kx  approaches  the  incident  material  wavenumber.  Note 
that  because  the  wavenumber  is  inversely  proportional  to  the  speed  of  sound,  the  trace 
sound  speed  on  the  surface  is  larger  than  the  speed  of  sound  in  air  for  all  non-grazing 
(' 9  <  90°)  angles. 


Figure  3.1.  Plane  wave  incidence  on  a  plane  surface.  An  incoming 
wave  with  k  =  uj/c  =  2n/X  at  an  angle  6  will  have  a  trace  wavenumber 
on  the  surface  of  kx  =  ksm(9). 


A  purely  real  value  of  k  is  typical  of  the  wavenumber  of  a  plane  wave,  thus  repre¬ 
senting  simple  phase  change  with  distance.  For  a  complex  k  value,  represented  as  the 
sum  of  the  real  and  imaginary  components  kr+jki,  the  plane  wave  equation  becomes: 

p  =  Aej$r+jni)-*  =  Aej%r-se-%i-s  ,  (3.4) 

with  the  imaginary  portion,  kt  =  A(fc),  used  to  represent  the  exponential  decay  or 
growth  of  an  evanescent  wave.  For  these  waves,  the  envelope  of  the  pressure  amplitude 
will  be  increasing  or  decreasing  (for  our  purposes  ”  decay /decrease”  will  henceforth 
be  used  in  lieu  of  accounting  for  both  possible  signs  of  fc*)  with  distance  (away  from  a 
defined  origin)  in  the  direction  of  propagation.  The  ratio  of  the  imaginary  component, 
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ki,  to  the  real  component,  kr,  defines  how  much  decay  occurs  over  each  period  of 
oscillation.  Over  a  wavelength  of  A  =  2n jkr  the  phase  increases  by  2i r;  when  the 
wavenumber  is  complex,  the  wave  decays  by  a  factor  of  e_27rfci/fcr  per  wavelength  (for 
a  1%  ratio  of  imaginary  to  real  component,  a  factor  of  0.93  per  wavelength). 

A  sample  of  a  wavetype  with  complex  wavenumbers,  as  viewed  in  the  direction  of 
propagation,  is  shown  in  Figure  3.2.  For  a  negative  value  of  ki,  the  wave  will  decay 
with  distance;  a  positive  value  of  ki  will  represent  an  increasing  wave  amplitude. 
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Figure  3.2.  The  effect  of  complex  wavenumbers  on  waveform  shape. 
The  undamped  wave  is  shown  in  the  left  plot.  A  1%  damping  (center 
plot)  will  cause  the  amplitude  to  decrease  over  distance,  while  a  1% 
negative  damping  (right  plot)  will  cause  the  amplitude  to  increase 
over  distance. 


3.1.2  Generation  of  Evanescent  Waves  Using  Complex  Angles 

A  complex  angle  of  propagation  can  be  used  to  represent  another  type  of  wave 
that  satisfies  the  wave  equation.  The  sine  or  cosine  of  this  complex  angle  will  result  in 
complex  wavenumber  components,  with  the  resulting  wavenumber  obeying  the  wave 
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equation.  The  complex  directional  wavenumbers  kx  and  kz,  are  determined  by  the 
sine  and  cosine  of  a  complex  angle  9  =  9r  +  j 9% : 


kx  =  k  sin(0r  ±  jOi )  =  k  sin  9r  cos (j9i)  ±  k  cos  9r  sin (j9i)  = 

k  sin  9r  cosh  9,  ±  jk  cos  9r  sinh  9i 
kz  =  k  cos (9r  ±  j9i )  =  k  cos  9r  cos (j9i)  =|=  k  sin  9r  sin (jOi)  = 

k  cos  9r  cosh  9i  =|=  jk  sin  9r  sinh  9Z  , 


(3.5) 


yielding  a  real  component  in  the  direction  of  propagation  and  an  imaginary  component 
that  is  normal  to  the  direction  of  propagation,  as  evidenced  by  the  cos (9r)  component 
in  the  imaginary  part  of  sine  of  the  complex  angle,  and  a  similar  sin(#r)  term  in  the 
cosine.  For  small  imaginary  angles,  the  hyperbolic  cosine  is  close  to  1;  therefore, 
the  resulting  directional  component  of  the  angle  (the  real  component  of  the  cosine 
of  the  complex  angle)  will  only  deviate  slightly  from  the  direction  of  propagation  for 
a  purely  real  angle.  The  imaginary  component  of  the  sine  or  cosine  of  the  complex 
angle  specifies  the  direction  of  decay  of  the  evanescent  wave,  with  the  hyperbolic  sine 
of  the  imaginary  component  indicating  the  rate  of  decay. 

Because  the  complex  angle  follows  the  standard  trigonometric  identity  sin2(0)  + 
cos2(0)  =  1,  the  use  of  the  complex  angle  allows  the  wave  propagating  at  that  angle 
to  continue  to  obey  the  wave  equation,  which  is  satisfied  by  the  resulting  directional 
vectors  being  equal  to  the  wavenumber  in  the  material: 


kl  +  k2z  =  k2  sin2(0)  +  k 2  cos 2 (9)  =  k2  ,  (3.6) 

which  is  an  essential  component  for  integrating  evanescent  waves  into  existing  theory. 
Note  that  a  complex  angle  acting  on  a  real  wavenumber  gives  a  very  different  result 
than  a  complex  wavenumber;  as  shown  in  Figure  3.3,  a  complex  angle  results  in  a  de¬ 
cay  perpendicular  to  the  direction  of  wave  propagation,  while  a  complex  wavenumber 
(and  a  real  angle)  decays  in  the  direction  of  propagation.  For  a  complex  value  of  the 
speed  of  sound,  c,  which  is  typically  used  to  model  damping  in  lossy  materials,  the 
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wave  equation  will  be  satisfied  by  a  complex  wavenumber  value.  When  using  com¬ 
plex  angles,  the  imaginary  component  of  the  wavenumber  transverse  to  the  direction 
of  propagation  represents  a  decay  per  wavelength  in  the  direction  perpendicular  to 
the  propagation  direction;  higher  wavenumber  magnitudes  will  result  in  faster  de¬ 
cay  per  distance  than  lower  wavenumbers  for  the  same  complex  angle  (note  that  the 
wavenumber  is  proportional  to  frequency). 
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Figure  3.3.  The  effects  of  a  complex  angle  of  propagation  (left)  com¬ 
pared  to  complex  wavenumber  (right).  The  left  pressure  distribution 
is  a  real  wavenumber  with  a  complex  angle,  causing  a  decay  normal  to 
the  direction  of  propagation.  The  right  distribution  shows  a  complex 
wavenumber  propagating  along  a  real  angle;  the  decay  is  along  the 
direction  of  propagation. 


When  incident  upon  a  surface,  an  evanescent  wave’s  trace  pressure  decays  with 
distance  over  the  surface.  Figure  3.4  compares  the  resulting  trace  pressure  distribu¬ 
tions  of  conventional  plane  and  evanescent  plane  wave  fields.  Since  our  concern  is 
with  interaction  occuring  at  a  material  interface,  generation  of  the  pressure  distribu¬ 
tion  through  a  volume  of  space  is  incidental;  in  Chapter  4  methods  for  generating 
these  distributions  on  the  surface  using  non-plane  wave  sources  will  be  explored.  The 
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use  of  complex  angles  is  a  convenient  way  to  express  how  these  pressure  distributions 
behave  through  the  interface. 


Figure  3.4.  The  resulting  trace  pressure  generated  by  incident  plane 
and  evanescent  waves  upon  a  surface.  The  upper  plots  show  propaga¬ 
tion  over  two-dimensional  space  above  the  incident  plane,  while  the 
bottom  plots  show  the  trace  pressure  distribution  at  the  z  =  0  surface. 
While  the  trace  pressure  of  the  30  degree  incident  plane  wave  shown 
on  the  left  does  not  decay  over  distance,  the  resulting  trace  pressure 
of  the  30  +  .1  j  degree  incidence  plot  on  the  right  decays  along  the 
plane  surface. 


3.1.3  Pressure  and  Velocity  Phase  Differences  in  Evanescent  Waves 

So  far  we  have  only  discussed  classification  of  evanescent  waves  by  their  change 
in  amplitude  over  distance.  However,  evanescent  waves  can  be  discerned  at  a  single 
point  in  space  by  the  difference  in  phase  between  the  wave’s  pressure  and  velocity. 
The  velocity  along  the  direction  of  propagation,  v,  of  a  plane  wave  is: 

_  1  dp 

jojp  dx 


(3.7) 
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Given  the  plane  wave  pressure  equation  in  Equation  3.1,  the  velocity  along  the  direc¬ 
tion  of  propagation  will  be: 


v  =  —p  ,  (3.8) 

up 

and  the  difference  in  phase  between  p  and  v  will  be  equivalent  to  the  phase  of  the 
wavenumber,  Zk.  For  a  complex  wavenumber,  k  =  kr+jki,  the  resulting  phase  angle, 
0,  is  equal  to: 


(j)  =  Zk  =  arctan 


(3.9) 


From  this  equation,  it  can  be  seen  that  there  is  no  phase  difference  between 
the  pressure  and  velocity  of  a  classical  plane  wave  with  k  =  kr,  and  a  90  degree 
difference  for  a  purely  imaginary  angle  of  propagation  k  =  ki.  For  a  complex  angle  of 
propagation,  the  phase  difference  between  the  velocity  and  pressure  of  a  evanescent 
wave  will  vary  between  0  and  90  degrees,  with  larger  phase  differences  corresponding 
to  a  more  signhcantly-decaying  wave. 

The  use  of  complex  angles  provides  a  way  to  simulate  waves  with  evanescent 
propagation  patterns,  and  can  be  used  to  express  refraction  characteristics,  which  is 
a  key  component  of  wave  transmission  across  an  interface. 


3.2  Evanescent  Wave  Transmission  Across  An  Interface 
3.2.1  Theory  and  Principles 

In  a  fluid,  which  cannot  sustain  shear  stress,  the  only  type  of  wave  that  can  propa¬ 
gate  is  a  longitudtinal  wave,  where  the  particle  velocity  of  the  wave  is  in  the  direction 
of  propagation.  Such  waves  are  also  known  as  compressional  waves  because  the  re¬ 
sulting  pressure,  p,  causes  compression  and  expansion  as  the  wave  propagates.  A  solid 
can  sustain  shear  stress  and  has  a  stress  vector,  <x,  instead  of  a  single  pressure;  there¬ 
fore,  an  additional  wavetype  with  speed  of  propagation,  6,  is  added  that  accounts  for 
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shear  motion  normal  to  the  direction  of  propagation.  In  a  three-dimensional  analysis, 
the  shear  wave  can  have  components  of  stress  and  velocity  in  independent  directions 
normal  to  the  direction  of  propagation;  for  our  purposes  in  a  two-dimensional  formu¬ 
lation,  the  shear  wave  stress  will  only  have  a  single  amplitude,  at-  The  compressional 
wave  will  have  a  stress,  cq,  which  is  equivalent  to  the  pressure  p  in  a  fluid.  Both 
compressional  and  shear  waves  are  present  in  a  solid  material,  and  both  must  be 
accounted  for  in  fluid-solid  propagation.  There  are  additional  wavetypes  that  may 
be  present  in  solid  materials,  such  as  Lamb  waves  that  cause  flexural  and  dilitational 
motion;  these  waves  result  from  the  coupling  of  the  pressure  and  shear  waves  in  con¬ 
strained  materials.  Because  these  are  geometry-dependent,  they  are  not  considered 
in  this  formulation. 

The  acoustical  impedance  of  a  material  is  the  ratio  of  stress  a  (for  fluids,  the 
scalar  —  p  may  be  used  to  designate  the  pressure;  note  that  a  is  conventionally  the 
negative  of  p)  over  the  particle  velocity,  v,  of  a  propagating  wave  in  the  material.  For 
a  solid  or  fluid,  the  characteristic  impedance  Z  is: 

Z  =  —  =  pc,  (3.10) 

v 

where  p  is  the  density  of  the  material  ( p0  and  c0  will  be  used  as  the  density  and 
speed  of  sound  of  air,  respectively,  with  the  ‘O’  subscript  representing  the  material 
properties  of  incident  acoustic  waves  in  air  throughout  this  document).  A  material 
that  supports  compressional  and  shear  waves  will  have  impedances  for  each  wavetype, 
proportional  to  their  propagation  speeds.  Table  3.1  shows  the  wave  speed  and  density 
for  several  common  materials. 

Ideal  energy  transmission  will  occur  between  materials  of  equivalent  impedance; 
for  interfaces  with  a  large  impedance  difference,  most  of  the  energy  in  a  wave  will 
reflect  back  into  the  incident  material.  For  typical  air-solid  interaction,  the  impedance 
difference  is  several  orders  of  magnitude,  and  a  very  small  amount  of  the  energy  is 
propagated  into  the  material  even  before  accounting  for  refraction  of  the  wave. 
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Table  3.1.  Density  and  wave  speed  of  materials  of  interest  for 
impedance  calculation. 


Material 

Density 

P 

kg/m3 

Compressional 

wave  speed 

c 

m/s 

Shear 

wave  speed 

b 

m/s 

Aii- 

1.2 

343 

N/A 

Water 

998 

1,480 

N/A 

Aluminum 

2,700 

6,110 

3,040 

Steel 

7,850 

5,850 

3,190 

Polycarbonate 

1,160 

1,700 

1,410 

Similar  to  a  material’s  wavenumber,  the  impedance  is  a  scalar  representing  the 
properties  in  the  material,  which  can  be  represented  in  a  propagating  wave  as  vector 
quantities  oriented  along  the  direction  of  propagation.  The  magnitude  of  a  wave’s 
impedance  vector  is  the  material  impedance,  and  each  component’s  amplitude  is  a 
function  of  the  angle  of  propagation.  In  a  fluid,  the  pressure  is  a  scalar  quantity,  and 
the  velocity  of  the  wave  is  oriented  along  a  vector;  therefore,  the  wave  impedance’s 
components  will  be  oriented  similarly  to  the  velocity’s.  For  a  surface  in  the  x  —  y 
plane  with  a  normal  vector  oriented  along  the  5  axis  (as  will  be  considered  in  the  rest 
of  this  formulation),  the  impedance  normal  to  the  material  is: 


z.=  '  = 


P 


|u  I  cos  6  cos  6 


(3.11) 


Impedance-matching  between  an  incident  wave  and  a  larger-impedance  material  can 
therefore  be  performed  by  changing  the  angle  of  the  incidence  of  the  wave.  However, 
for  most  air-solid  interfaces  with  large  material  impedance  differences,  the  angle  that 
will  match  the  normal  impedances  is  extremely  close  to  grazing,  making  it  physically 
impractical  to  implement.  Even  with  the  ideal  angle  of  incidence  required  to  obtain 
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an  impedance-matched  system,  the  angle  of  propagation  for  the  transmitted  wave 
due  to  refraction  at  the  interface  can  reduce  the  effective  energy  propagation  into  the 
material. 

A  difference  in  the  speed  of  sound  between  two  materials  causes  refraction  across 
the  interface.  The  speed  of  any  propagating  waves  will  be  constant  parallel  to  the 
surface:  i.e.,  the  trace  wave  speed  is  common  between  the  incident  and  transmitted 
media.  Therefore,  for  a  wave  incident  at  an  angle  to  the  surface:  the  resulting  angles 
of  the  compressional  and  shear  waves  can  be  calculated  from  Snell’s  law: 

Co  _  ci  _  bi 
sin(#0)  sin(di)  sin(7i) 

where  Ci  and  0\  are  the  compressional  wave  speed  and  propagation  angle  of  the 

transmitted  wave,  respectively,  and  by  and  71  are  their  shear  equivalents.  The  effect 

of  refraction  across  the  surface  is  that  the  angle  of  propagation  into  a  higher-speed 

material  will  be  a  larger  than  that  of  the  incident  wave. 

For  a  wave  incident  upon  a  plane  surface,  the  angle  of  incidence  determines  the 
kx  component  on  the  surface;  this  component  is  the  same  for  both  the  incident  and 
transmitted  waves,  so  kx  =  k\x  =  kox :  i.e., 

k0  sin  9(j  =  kisinOi  =  Ky  sin  7^  .  (3.13) 

The  resulting  k\z  of  the  transmitted  wave  is  determined  by  calculating  the  re¬ 
maining  component  of  the  wavenumber  vector: 

kyz  =  \jk\  —  k2x  =  k\  cos  0\  ,  (3-14) 

and  the  reduction  in  kz  due  to  the  wavenumber  changes  across  an  interface  can  be 
seen  as  another  aspect  of  refraction. 

The  formulation  of  kz  as  the  remaining  component  after  removing  the  trace 
wavenumber  is  helpful  to  understand  subsonic  propagation.  When  the  trace  wavenum¬ 
ber  is  higher  than  the  material  wavenumber,  the  trace  speed  of  sound  on  the  surface 


(3.12) 
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is  lower  than  the  material  wave  speed.  In  this  case,  the  wave  component  propagating 
into  the  material  is  supersonic  (faster  than  the  material  wave  speed),  and  will  prop¬ 
agate  like  a  conventional  plane  wave.  If  the  trace  wavenumber  exceeds  the  material 
wavenumber,  the  resulting  wave  into  the  material  is  subsonic;  the  wave  will  propa¬ 
gate  parallel  to  the  interface,  but  will  be  evanescent  and  create  a  near-held  pressure 
distribution  along  the  interface.  A  subsonic  wave  producing  an  evanescent  pressure 
distribution  will  not  propagate  energy  into  the  material.  In  the  supersonic  region,  the 
angle  of  propagation  is  entirely  real;  the  shift  from  supersonic  to  subsonic  propagation 
is  the  shift  from  a  real  to  complex  angle  (where  the  real  component  of  the  angle  will 
always  equal  90°).  The  critical  angle,  9oc,  for  an  incident  wave  where  the  resulting 
wave  into  the  material  becomes  subsonic  is: 

9qc  =  arcsin  —  .  (3.15) 

Cl 

For  air-solid  interfaces,  this  angle  is  very  small;  for  example,  for  a  ratio  of  wave 
speeds  ci/c0  =  10,  the  critical  angle  is  5.7  degrees.  For  angles  greater  than  the 
critical  angle,  the  transmitted  angle  into  the  material  is  complex,  and  the  resulting 
normal  wavenumber  into  the  transmitted  material  will  be  imaginary. 

The  evanescent  waves  generated  by  complex  angles  increase  energy  propagation 
by  decreasing  the  refracted  real  angle  and  allowing  for  wave  propagation  at  angles 
higher  than  the  critical  angle.  Figure  3.5  shows  the  change  in  both  real  and  complex 
components  of  a  refracted  angle,  9 1,  from  a  complex  incident  angle,  9q,  across  an 
interface  with  C\/cq  =  10.  The  addition  of  an  imaginary  component  to  the  angle 
eliminates  the  critical  angle  criterion.  The  refraction  of  a  complex  angle  will  asymp¬ 
tote  towards  90°;  adding  only  a  small  imaginary  component  to  the  complex  angle 
causes  significant  motion  away  from  grazing.  For  example,  an  A(di)  value  of  1  causes 
the  refracted  angle  at  30°  to  shift  from  grazing  to  37°.  However,  the  imaginary  com¬ 
ponent  increases  with  refraction  as  well;  the  addition  of  the  imaginary  component  to 
the  angle  causes  an  increase  in  the  imaginary  component  of  the  propagating  angle 
from  1.3  (which  occured  at  grazing)  to  1.8. 
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Imaginary  Component  of  6^  (radians)  Imaginary  Component  of  0Q  (radians) 

Figure  3.5.  Complex  angle  refraction  of  the  transmitted  wave 
9\  =  6r  —  jOi  across  an  interface  with  ci/co  =  10.  The  left  plot  is 
K(arcsin(ci/c0  sin  0O)),  while  the  right  plot  is  3?(ci/c0  arcsin(sin  90)). 


Even  with  a  reduction  in  the  refracted  angle  due  to  the  use  of  complex  incident 
angles,  the  increase  of  imaginary  components  will  create  additional  decay  in  the  prop¬ 
agating  wave.  By  calculating  the  energy  decay  with  distance  across  the  interface,  we 
can  more  accurately  see  the  effects  of  the  evanescent  wave.  The  calculation  of  in¬ 
tensity  would  include  both  the  pressure  transmission  coefficient  and  the  effect  of  the 
complex  angle  in  wave  transmission. 

The  calculated  intensity  transmission  coefficients  are  functions  of  the  material 
impedances  and  incoming  angles  of  the  wavetypes  in  the  system.  Considering  both  of 
these  factors  will  yield  a  full  accounting  of  the  effects  of  evanescent  waves  on  energy 
transmission  across  an  interface. 
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3.2.2  Fluid-Fluid  Interface 

We  will  begin  by  considering  a  fluid-fluid  interface,  since  it  is  a  simpler  system  with 
only  one  wavetype  in  each  material.  Our  example  system  will  be  an  air-fluid  interface 
where  the  theoretical  fluid  has  speed  of  sound  c\  =  10co  and  density  p\  =  lOOOpo- 
At  the  fluid-fluid  interface  there  are  three  waves  to  consider:  the  incident  wave 
incoming  at  an  angle  90  of  unit  amplitude,  the  reflected  wave  in  the  incident  fluid 
at  an  angle  —  90  with  amplitude  R ,  and  the  transmitted  wave  at  angle  9\,  computed 
from  Snell’s  law  as  a  function  of  the  incoming  angle  and  two  fluids’  speed  of  sound, 
and  amplitude  T.  For  an  interface  in  the  x  —  y  plane  with  a  normal  z,  the  material 
states  across  the  interface  remain  constant: 

Pi(0)  =  p0(0)  ;  uu(0)  =  v0z  .  (3.16) 

Considering  the  incident,  reflected,  and  transmitted  waves  in  these  formulations 
yields: 


T  —  1  +  R;  - I-^-  =  — - (L-  ,  (3.17) 

Z\l  cos(6»i)  Z0/cos(6»0) 

and  these  two  boundary  conditions  can  be  used  to  solve  for  the  transmitted  wave 
amplitude  T : 


‘IZxj  cos(^i) 


T  = - (3.18) 

Z1/ 008(01)  +  V  cos(0o) 

Note  that  each  additive  term  in  the  equation  is  equivalent  to  the  normal  impedance  in 
the  fluids.  The  amplitude  of  T  with  incident  angle  in  a  classical  system,  as  shown  in 
Figure  3.6,  is  nearly  2  for  all  incident  angles  due  to  the  dominance  of  the  impedance 
of  water.  Above  the  critical  angle,  T  is  equal  to  2,  but  the  transmitted  angle  will  be 
imaginary  and  no  energy  will  propagate.  At  incident  angles  extremely  close  to  graz¬ 
ing,  the  influence  of  the  incident  impedance  will  reduce  the  transmission  coefficient; 
however,  there  is  still  no  energy  propagation  into  the  material. 
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Figure  3.6.  Pressure  transmission  coefficient  T  for  plane  wave  prop¬ 
agation  into  a  fluid. 


For  complex  90,  the  pressure  transmission  coefficient  will  decrease  with  increasing 
imaginary  angle  component.  A  plot  of  the  transmitted  pressure  amplitude  is  shown 
in  Figure  3.7,  which  shows  that  the  maximum  values  of  transmission  occur  at  real, 
grazing  angles. 

The  pressure  transmission  coefficient  does  not  account  for  the  resulting  velocity 
of  the  transmitted  wave,  which  is  necessary  to  discern  the  energy  propagating  into 
the  system.  Development  of  an  intensity  transmission  coefficient,  which  will  properly 
account  for  the  energy  transmitting  across  the  interface  due  to  both  the  pressure 
transmission  coefficient  and  the  complex  angle  of  propagation,  is  therefore  necessary 
to  compute  the  sound  power  transmitted. 

The  intensity  of  a  wave  is  computed  from  the  pressure  and  velocity  at  that  point: 

I  =  ^3?(p(T)u*(T))  , 


(3.19) 
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Figure  3.7.  Pressure  transmission  coefficient  T  for  fluid-fluid  system. 


where  v*  is  the  complex  conjugate  of  the  velocity.  This  relationship  accounts  for  both 
the  pressure  and  velocity  amplitudes  (which  are  related  by  the  impedance)  as  well  as 
the  phase  relationship  between  the  pressure  and  velocity. 

The  velocity  in  the  z  direction  of  a  plane  wave  is  given  by  Equation  3.7;  for  a  plane 
wave  with  the  pressure  distribution  shown  in  Equation  3.1  (note  that  the  wavenumber 
vector  has  been  expanded  to  it’s  component  terms),  the  resulting  velocity  is: 

vz  =  —p  ,  (3.20) 

ujp 

and  the  resulting  general  equation  for  intensity  in  the  plane  wave  normal  to  the  surface 
is: 


h  =  -m  — 


(kz  A 

— PP  ■ 
\up  J 


(3.21) 


For  a  real  incident  angle,  the  intensity  equation  simplifies  to  the  ratio  of  the  mean 
square  pressure  over  the  material  impedance: 
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Iz  =  ^-c°s6  ■  6  =  ^(6)  ,  (3.22) 

where  \A\2  is  the  magnitude  squared  of  the  complex  amplitude,  \A\2  =  AA*.  For 
subsonic  wave  propagation  (a  complex  angle  with  a  90  degree  real  component),  the 
wavenumber,  kz.  in  the  z  direction  will  be  purely  imaginary,  and  the  intensity  in  the 
wave  will  be  equal  to  zero.  For  a  complex  angle  of  incidence,  9o  =  90r  +  j0Oi,  the 
resulting  intensity  in  the  wave  will  be: 


L  = 


_  c2Q(kzz-2kxx ) 

2up 


(3.23) 


The  intensity  transmission  coefficient  across  an  interface  in  the  z  direction,  It, 
resulting  from  a  unit  incident  pressure  amplitude  is: 


=  w  =  pm.sv,-i...i 

/o(*o)  ^{k,QZ)pi 

and  for  the  intensity  transmission  directly  across  the  interface,  Zo  =  Z\ 
intensity  transmission  coefficient  can  be  simplihed  to: 


(3.24) 
0,  the 


/i(0)  \T\2^klz)p0 

T  /o(0)  %(k0z)pi 

with  the  transmission  coefficient  TC  in  decibels  equal  to: 


(3.25) 


TC  =  10  log(Fr)  ,  (3.26) 

which  is  the  inverse  of  the  transmission  loss. 

The  inclusion  of  an  imaginary  angle  component  become  beneficial  when  consider¬ 
ing  intensity  transmission  across  the  interface.  The  intensity  transmission  coefficient 
as  a  function  of  real  and  complex  angle  is  plotted  in  Figure  3.8.  The  plot  shows 
increased  intensity  transmission  with  increasing  imaginary  angle  component;  for  ex¬ 
ample,  an  increase  in  imaginary  angle  component  from  0.001  to  0.3  causes  23.6  dB  of 
increase  to  the  transmission  coefficient  at  a  30  degree  incident  angle. 
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Figure  3.8.  Intensity  transmission  coefficient  TC  variation  for  a  fluid- 
fluid  system  with  a  complex  angle  of  incidence. 

3.2.3  Fluid-Solid  Interface 

Wave  propagation  at  a  fluid-solid  interface  will  generate  a  reflected  pressure  wave 
in  the  incident  fluid  and  transmitted  shear  and  pressure  waves  in  the  transmitted 
solid,  with  the  shear  amplitude  B  and  a  wavenumber  k.  At  the  fluid-solid  material 
interface,  pressure  and  velocity  are  conserved  across  the  boundary.  Because  a  fluid 
wave  cannot  sustain  shear,  the  shear  stress  on  the  surface  is  zero.  Note  that  this 
does  not  mean  the  shear  wave  amplitude  is  zero;  the  compressional  wave  will  induce 
a  shear  stress  on  the  surface  when  propagating  at  a  non-normal  angle  and  the  shear 
wave  stress  is  generated  to  cancel  out  this  surface  stress.  The  three  state  relationships 
at  the  interface  are: 


Po(0)  =  -az( 0)  ,  0  =  o-x(0)  ,  uO2(0)  =  vg(0)  . 


(3.27) 


For  an  incident  wave  with  unit  pressure  amplitude,  the  resulting  transmitted  pressure 
amplitudes  into  the  material  (as  calculated  by  Brekhoshovikh  [3])  is  proportional  to 
the  angles  of  propagation  and  the  normal  impedances  at  the  surface: 
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7  -  pC  ■  7  -  PlCl  ■  7  -  Plbl 
cos(90 )  ’  lz  cos(ffi)  ’  lx  cos(7i) 

For  a  unit  amplitude  incident  wave,  the  amplitude  of  the  normal  stress  propagating 

into  the  material  is  the  transmission  coefficient  Tg\ 


(271)  +  Zz 

while  the  resulting  shear  wave  transmission  coefficient  Tg  is: 

T  _ _ 2Zlx  sin(27i) _ 

5  Zlz  cos2  (271)  +  Zlx  sin2  (271)  +  Z2  ' 

The  pressure  transmission  coefficient  for  a  real  angle  has  several  distinct  charac¬ 
teristics  due  to  the  interplay  of  the  compressional  and  shear  waves.  Figure  3.9  shows 
the  pressure  transmission  coefficient  in  logarithmic  form  ( TC  =  20  log (T^g)  for  an 
interface  with  c\/cq  =  10,  61/co  =  7  and  density  ratio  pi/po  =  1000.  There  are  sev¬ 
eral  peaks  in  the  values  that  would  appear  to  indicate  avenues  of  energy  transmission. 
The  small  spike  in  Tg  at  5.7  degrees  is  due  to  the  critical  angle  of  the  compressional 
wave, which  causes  a  maxima  in  the  numerator  of  Tg.  The  dip  at  5.8  degrees  is  caused 
by  a  minima  in  the  numerator  when  cos(27i)  is  equal  to  zero.  The  critical  angle  for 
the  shear  wave  is  8.2  degrees,  where  there  is  a  discontinuity  in  both  the  Tg  and  Tg 
value.  The  peak  near  9.4  degrees  is  due  to  impedance  cancellation  (which  occurs  af¬ 
ter  the  critical  angle  and  71  is  complex)  between  the  transmitted  pressure  and  shear 
wave,  minimizing  the  denominator  common  to  both  terms.  However,  most  of  these 
features  are  at  angles  beyond  the  critical  angles  for  both  the  compressional  and  shear 
waves,  and  there  is  very  small  velocity  into  the  material. 

Intensity  in  the  solid  is  calculated  in  a  way  similar  to  Equation  3.23,  with  the 
appropriate  wave  amplitude  and  wavenumber  in  place  for  the  compressional  or  shear 
wave: 


(3.30) 


Tr  = 


2 Zlz  cos (271 


Zlz  COS2  (271)  +  Zlx  sin2 


(3.29) 


Il,z  = 


\A\2^{k 


z)  2^s(kzz—2kxx) 


2up 


(3.31) 
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Figure  3.9.  Transmission  coefficients,  201og(Ti)  and  201og(7s),  for 
a  system  with  wave  speed  ratio  ci/co  =  10  and  density  ratio  pi/po  = 


1000. 


Is,z  ~ 


_  2^{kzz-2k, 


x) 


2c op 


(3.32) 


Because  intensity  is  a  measure  of  energy,  it  should  balance  across  the  interface: 
i.e. ,  the  total  power  incident  on  the  surface  minus  the  intensity  reflected  above  the 
surface  should  equal  the  intensity  transmitted  below  the  surface  in  longitudinal  and 
shear  waves.  The  model  above  can  be  verified  by  accounting  for  this  intensity.  The 
magnitude  of  the  reflected  wave  component,  R ,  is  calculated  as: 


R  _  Z\z  cos2(27i)  +  Zlx  sin2  (271)  - 

Ziz  cos2  (271)  +  Zlx  sin2  (27! )  +  Zz  ' 

The  reflected  intensity  travels  away  from  the  interface;  therefore,  the  net  intensity 
above  the  interface  will  be  the  intensity  of  the  incident  wave  minus  the  intensity  of 
the  reflected  wave.  Figure  3.10  shows  the  total  intensity  above  and  below  the  fluid 
interface  for  both  a  classical  and  evanescent  wave.  That  the  intensity  is  the  same 
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above  and  below  the  interface  across  the  incident  angle  range  serves  as  verification 
that  the  model  is  properly  accounting  for  energy  in  the  system. 


Figure  3.10.  Net  normal  intensity  above  and  below  the  fluid-solid 
interface,  with  the  incident  and  reflected  waves  contributing  to  the 
intensity  above  the  interface  and  the  longitudinal  and  shear  transmit¬ 
ted  waves  contributing  to  energy  below  the  interface. 


The  intensity  transmission  coefficient  across  the  material  interface  used  previously 
(ci/co  =  10,  bi/co  =  7  and  pi/po  =  1000)  is  shown  below  in  Figure  3.11.  This  inten¬ 
sity  transmission  is  that  immediately  across  the  interface,  and  will  not  represent  the 
effect  of  the  evanescent  component  penetrating  further  into  the  material.  The  lowest 
complex  incidence  angle  value  shown  is  0.001,  as  a  purely  real  value  will  be  equal  to 
zero  beyond  the  critical  angle.  The  use  of  complex  angles  expands  the  transmission 
peaks  (such  as  those  at  10  Hz)  while  allowing  for  energy  transmission  across  the  entire 
spectrum.  While  the  complex  angle  does  not  provide  a  substantial  benefit  below  the 
critical  frequency  (the  warm  area  to  the  far  left  of  the  plot),  expansion  of  the  realm 
of  intensity  transmission  beyond  the  near-normal  angles  will  allow  for  more  varied 
sources  of  excitation. 
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Figure  3.11.  Intensity  transmission  coefficients  at  1000  Hertz  across 
the  interface  for  the  compressional  (left)  and  shear  (right)  waves  due 
to  a  complex  angle  of  the  incident  wave. 


The  intensity  transmission  coefficient  calculated  at  a  depth  of  one  meter  below 
the  material  interface  is  shown  in  Figure  3.12.  This  transmission  coefficient  shows  the 
negative  effects  associated  with  the  imaginary  angle  of  refraction.  There  is  significant 
intensity  lost  into  the  system,  particularly  above  the  critical  angle,  and  options  for 
ideal  intensity  transmission  are  limited  to  smaller  angles.  However,  there  is  still 
transmission  above  the  critical  angle,  and  the  pressure  transmission  peak  near  10 
degrees  provides  significant  transmission. 

The  use  of  complex  angles  to  generate  evanescent  waves  has  been  shown,  in  theory, 
to  increase  energy  transmission  across  the  air-solid  interface.  While  using  complex 
angles  introduces  additional  imaginary  components  to  the  refracted  angle  of  the  prop¬ 
agating  wave,  the  change  in  real  angle  of  propagation  will  cause  energy  to  propagate 
further  into  the  material.  An  approach  for  intensity  propagation  taking  into  account 
the  pressure  transmission  and  the  resulting  intensity  due  to  the  angle  of  transmission 
must  be  formulated  to  allow  for  a  means  to  see  the  effect  of  complex  angles  on  energy 
transmission  across  the  interface. 
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Figure  3.12.  Intensity  transmission  coefficients  for  a  depth  of  1  me¬ 
ter  into  the  material  for  the  compressional  and  shear  waves  due  to 
complex  angles. 


The  use  of  complex  angles  can  increase  intensity  transmission  at  a  high-impedance 
interface,  particularly  beyond  the  critical  angle  of  sound  transmission.  By  using  the 
evanescent  components  of  a  monopole,  we  can  show  that  these  theoretical  considera¬ 
tions  can  be  seen  in  experimental  verification. 

3.3  The  Monopole  as  an  Evanescent  Wave  Source 

Evanescent  pressure  distributions  with  well-defined  characteristics  are  not  gener¬ 
ated  by  typical  sound  sources.  While  some  sources  may  generate  evanescent  waves, 
they  will  be  unique  to  the  source  and  difficult  to  control.  Therefore,  it  is  helpful  to  ex¬ 
plore  evanescent  wave  properties  of  a  simple  source,  the  acoustic  monopole,  to  provide 
a  basis  for  experimental  exploration  of  the  theory  of  evanescent  wave  transmission. 
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3.3.1  Theory 

An  acoustic  monopole  is  a  fundamental  source,  in  which  pressure  generation  from 
a  single  point  radiates  equally  in  all  directions.  The  pressure  at  a  distance  R  from 
the  source  is  equal  to: 


pjk0R 

P  =  ,  (3-34) 

which  shows  a  1  /R  decay  in  pressure  over  distance.  Figure  3.13  shows  the  pressure 
of  a  monopole  in  space  and  along  the  radius  from  the  monopole. 


Figure  3.13.  Monopole  pressure  distribution  in  the  plane  containing 
the  monopole,  showing  the  real  component  of  the  pressure  for  a  1000 
Hz  monopole  in  air  (a  wavenumber  of  18.3  radians  per  meter),  in 
two-dimensional  (left)  and  in  one- dimensional  space  (right)  along  the 
radius.  The  maximum  pressure  peak  at  r  =  0  has  been  cut  off  at  10 
Pa  to  show  detail. 


When  the  radiated  pressure  distribution  is  projected  onto  a  plane  offset  at  a 
distance  Zq,  the  resulting  distance  from  the  origin  R  will  be: 


R  = 


X2  +  y2  +  Zg  , 


(3.35) 
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and  the  radius  will  increase  from  a  minimum  value  of  Zq  along  the  plane.  A  larger  z0 
value  will  cause  a  reduced  increase  in  R  per  length  along  the  plane,  which  will  reduce 
the  amount  of  decay  per  distance  in  the  trace  pressure  distribution.  Furthermore,  the 
trace  wavelength  along  the  plane  will  decrease  to  the  material  wavelength  2n/k  as 
R  increases;  a  higher  zq  value  will  cause  less  variation  in  the  trace  wavelength.  The 
result  of  the  decrease  in  decay  rate  and  trace  wavelength  will  cause  the  monopole 
to  approximate  a  plane  wave.  Sample  trace  pressure  distributions  for  both  a  small 
and  large  standoff  distance  are  shown  in  Figure  3.14.  The  10  centimeter  standoff  has 
significant  decay  near  the  origin  and  shorter  wavelengths,  while  the  2  meter  standoff 
has  less  decay  and  the  initial  wavelengths  are  closer  to  the  material  wavelength. 


Figure  3.14.  Trace  pressure  distributions  (real  component  of  pres¬ 
sure)  of  a  monopole  at  a  distance  of  10  and  200  centimeters  for  a 
wavenumber  of  18.3  radians  per  meter.  The  shorter  standoff  distance 
displays  similar  properties  to  an  evanescent  wave. 


Due  to  the  1/R  decay  (compared  to  an  exponential  decay)  and  the  change  in 
wavelength,  a  monopole  pressure  distribution  is  not  equivalent  to  an  evanescent  plane 
wave.  However,  it  does  display  evanescent  properties. 
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The  wavenumber  spectrum,  which  can  be  calculated  by  the  Fourier  transform  of 
the  spatial  pressure  distribution: 


P{k)  =  F{p{x))i,  (3.36) 

of  a  monopole’s  trace  pressure  distribution  consists  of  distinct  regions  of  evanescent 
and  non-evanescent  waves.  The  energy  in  the  wavenumber  region  below  ko  is  super¬ 
sonic,  while  the  energy  above  ko  is  subsonic  and  generates  evanescent  waves.  Figure 
3.15  shows  the  wavenumber  spectrum  for  the  two  pressure  distributions  shown  in 
Figure  3.14.  With  increasing  standoff  distance,  the  resulting  WN  spectra  has  less 
energy  in  the  evanescent  region.  Figure  3.16  shows  Q ,  the  ratio  of  evanescent  com¬ 
ponent  to  non-evanescent  component  in  the  trace  pressure  distribution  as  a  function 
of  normalized  standoff  distance: 


/~P2(Mo) 


Qikzo)  =  ^ - — —  (3.37) 

I0°P2(k,z0) 

A  monopole  produces  significant  evanescent  components,  and  these  components 
can  be  decreased  by  increasing  the  standoff  distance  between  the  monopole  and  the 
incident  plane.  Using  these  properties  of  the  monopole,  the  effects  of  evanescent 
waves  on  transmission  of  energy  into  the  plate  can  be  seen. 


3.4  Experimental  Measurement  of  Monopole  Evanescent  Components 

The  monopole  provides  a  source  of  evanescent  waves,  which  may  propagate  ad¬ 
ditional  energy  into  higher  impedance  materials.  The  pressure  distribution  of  a  sim¬ 
ulated  monopole  source  was  inspected  for  the  presence  of  evanescent  waves  using 
acoustical  holography.  Once  the  evanescent  waves  were  verified,  the  distribution  was 
used  incident  upon  a  panel  to  see  the  effects  of  the  evanescent  wave  transmission  into 
the  solid  material. 
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Figure  3.15.  Wavenumber  spectra  of  the  trace  pressure  distribution 
from  a  monopole  radiating  at  100  Hz  for  standoff  distances  of  10  and 
200  centimeters.  The  peak  in  the  plot  is  at  k0 ;  energy  above  this  peak 
is  subsonic,  with  the  closer  standoff  having  more  subsonic  energy  than 
the  further  standoff. 

3.4.1  Test  Setup  and  Methodology 

The  panel  used  for  testing  was  a  9.4  millimeter  (nominally  0.375  inch)  thick  poly¬ 
carbonate  panel.  The  panel  was  60  centimeters  long  by  50  centimeters  high,  with 
holes  drilled  in  each  corner  at  12.7  millimeters  (0.5  inch)  from  the  edges.  These 
holes  were  used  to  suspend  the  plate  in  a  113  centimeter  square  wooden  frame.  This 
mounting  method  was  used  to  minimize  the  effects  of  mounting  on  the  vibration  of 
the  plate,  as  well  as  to  minimize  reflections  from  the  mounting.  The  plate  was  sanded 
on  one  side  to  produce  an  opaque  surface  for  improved  laser  vibrometer  focusing  and 
measurement. 

Vibration  properties  for  the  plate  were  measured  by  exciting  the  panel  with  a 
PCB  712A02  stack  actuator  equipped  with  a  100  gram  add-on  mass.  The  resulting 
force  into  the  system  was  measured  with  a  PCB  208A02  force  transducer  mounted 
between  the  actuator  and  the  panel  surface,  and  was  recorded  using  the  data  acqui- 
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Figure  3.16.  Ratio  Q  of  evanescent  to  subsonic  to  supersonic  com¬ 
ponents  in  the  trace  pressure  distribution  of  a  monopole  incident  on 
a  plane  as  a  function  of  normalized  monopole  standoff  distance.  The 
increase  in  standoff  causes  a  reduction  in  evanescent  energy. 


sition  system  from  the  laser  vibrometer.  This  actuactor-transducer  arrangement  was 
mounted  in  the  corner  of  the  panel  during  vibration  measurements  to  classify  modes 
and  wavetypes  in  the  panel. 

Vibration  was  measured  using  the  Polytec  PSV  400  scanning  laser  vibrometer  and 
data  acquisition  system.  The  vibrometer  allows  for  automated  scanning  of  points 
along  a  surface.  Vibration  was  recorded  from  0  to  4  kilohertz,  along  with  the  output 
from  either  the  force  transducer  or  microphone  near  the  source  (depending  on  the 
test). 

The  vibrometer  data  acquisition,  which  was  used  for  all  vibration  measurements, 
measured  a  bandwidth  of  4,000  Hertz  at  2.5  Hz  resolution,  with  50  averages  per  point. 
A  total  of  221  points  were  measured,  constituting  13  rows  and  17  columns  across  a 
54  by  42  centimeter  rectangle  on  the  panel.  The  resulting  measurement  parameters 
allow  for  a  wavenumber  resolution  of  11.63  radians  per  meter  (rad/m)  from  -69.81 
rad/m  to  58.18  rad/m  in  the  horizontal  direction  and  a  resolution  of  15.14  rad/m 
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from  -121.12  to  105.98  rad/m  in  the  vertical  direction.  More  information  about 
wavenumber  decomposition  and  analysis  can  be  found  in  Section  5.2. 

The  monopole  pressure  distribution  was  generated  using  the  Herrick  Laboratories 
half-monopole  source.  This  source  is  a  1  meter  long,  2.54  cm  (1  inch)  diameter  copper- 
lined  PVC  pipe  attached  to  a  CIE  PD30T  30W  compression  driver  with  a  2.54  cm 
(1  inch)  speaker  diameter.  The  length  of  PVC  pipe  ensures  plane  wave  propagation 
in  the  tube,  that,  when  radiating  from  the  open  end  of  the  pipe,  produces  a  near¬ 
monopole-like  pressure  distribution  in  the  hemispherical  space  centered  at  the  pipe 
opening.  The  validity  of  the  resulting  pressure  distribution  will  be  discussed  in  the 
following  section.  While  this  source  will  produce  significant  harmonics  at  the  half¬ 
wavelength  frequencies  of  the  pipe,  it  is  suitable  for  single-frequency  analysis.  The 
source  was  driven  using  a  General  Radio  1381  Random-Noise  Generator  at  5  kHz  and 
a  QSC  1100  amplifier. 

Measurement  of  sound  radiation  from  the  half-monopole  source  was  performed 
using  a  Briiel  and  Kjaer  4189  microphone  placed  immediately  behind  the  pipe  opening, 
approximately  2  centimeters  from  the  surface.  While  this  location  is  less-than-ideal 
for  measuring  the  output  spectra  of  the  monopole,  it  was  chosen  so  as  not  interfere 
with  the  radiation  pattern.  Additional  testing  showed  that  the  spectra  obtained  from 
this  location  matched  closely  with  the  spectra  obtained  from  the  pipe  opening.  The 
microphone  was  connected  to  a  Briiel  and  Kjaer  2169  preamp  and  5935  power  supply, 
and  into  the  data  acquisition  system  for  the  vibrometer. 

For  classification  of  the  monopole’s  sound  radiation,  the  pressure  distribution 
was  measured  with  a  64  microphone  holography  array,  consisting  of  TMS  T130C21 
ICP  microphones  in  an  8  by  8  square  pattern,  with  each  microphone  separated  from 
its  neighbors  by  10  centimeters.  The  array  provides  wavenumber  resolution  of  8.97 
radians  per  meter  from  -35.90  rad/m  to  26.93  rad/m.  In  air,  35.9  radians  per  meter 
occurs  at  1,960  Hertz;  the  array  should  provide  sufficient  range  to  properly  visualize 
wavenumber  spectra  at  frequencies  below  that.  The  array  output  was  recorded  using 
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a  VXI  data  acquisition  system  with  custom  MATLAB  software  for  measurement  and 
analysis. 

Testing  was  performed  at  distances  between  35  and  145  millimeters  from  micro¬ 
phone  number  43,  located  in  the  6th  row  (from  the  top)  and  3rd  column  (from  the 
left)  of  the  array.  The  source  was  measured  off-center  to  allow  for  additional  distance 
to  measure  the  decay  in  the  sound  field;  a  centered  microphone  would  not  have  the 
same  range  of  distance  data.  A  schematic  of  the  setup  is  shown  in  Figure  3.17. 
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Figure  3.17.  Setup  and  geometry  for  the  monopole  pressure  radiation 
using  the  holography  array. 


A  schematic  of  the  setup  for  measuring  acoustical  excitation  of  the  panel  is  shown 
in  Figure  3.18,  and  a  photo  of  the  experimental  setup  is  shown  in  Figure  3.19.  The 
monopole  was  positioned  at  the  center  of  the  panel  and  tested  at  various  distances 
from  2.2  centimeters  to  20.2  centimeters.  The  resulting  vibration  was  recorded  us¬ 
ing  the  vibrometer,  and  the  resulting  sound  propagation  was  measured  using  the 
microphones. 
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Figure  3.18.  Setup  for  panel  excitation  from  a  monopole.  The 
monopole  produces  a  pressure  distribution  on  the  panel,  and  the  re¬ 
sulting  vibration  is  measured  using  the  laser  vibrometer. 


3.4.2  Results 

For  the  purposes  of  brevity  and  clarity,  all  results  will  be  presented  at  1,000  Hertz. 
This  frequency  is  appropriate  because  it  does  not  correspond  to  a  mode  of  vibration 
in  the  plate,  and  is  representative  of  many  other  regions  of  plate  vibration. 

Two  ’’slices”  of  the  wavenumber- frequency  data,  at  the  center  of  the  horizontal 
(left)  and  vertical  (right)  wavenumber  domains,  are  shown  in  Figure  3.20.  The  vibra¬ 
tion  of  the  plate  is  dominated  by  a  dispersive  wave  with  a  wave  speed  of  approximately 
200  m/s  at  1000  Hz;  this  wave  will  be  subsonic  in  air.  The  wavenumber  spectrum 
of  vibration  is  shown  in  Figure  3.21,  showing  the  concentration  of  energy  in  an  oval 
with  major  axes  at  45.42  and  34.91  radians/meter  (although  the  precise  locations  of 
these  peaks  are  within  the  resolution  of  the  wavenumber  measurement). 
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Figure  3.19.  Experimental  setup  for  monopole-excited  panel.  The 
monopole  source  is  the  forefront  object,  with  the  driver  (light  blue) 
closest  to  the  camera.  The  panel  is  mounted  in  the  wooden  frame;  the 
white  surrounding  material  is  Thinsnlate,  which  is  wrapped  around 
the  frame  to  reduce  reflections  from  the  frame.  The  experiment  took 
place  in  the  hcmianechoic  chamber,  which  reduces  reflections  from 
other  sources. 


The  pressure  spectrum  of  the  monopole  shows  a  1  / R  decay  with  distance  (from 
the  monopole  location),  as  seen  by  comparing  the  pressure  at  each  measured  point  as 
a  function  of  the  radius  in  Figure  3.22.  The  wavenumber  spectrum  of  the  monopole 
pressure,  shown  in  Figure  3.23  shows  energy  in  a  circle  with  radius  of  approximately 
18  radians  per  meter,  coinciding  with  theory.  In  addition,  the  pressure  and  intensity 
from  the  monopole  decay  with  the  radius  (from  the  monopole)  and  radius  squared, 
respectively. 

The  ratio  of  evanescent  component  to  non-evanescent  component  (similar  to  what 
is  shown  in  Figure  3.16  for  Equation  3.37)  for  the  half-monopole  source  is  shown  in 
Figure  3.24.  Similar  to  a  theoretical  monopole,  the  half-monopole  source  shows  a 
reduction  in  evanescent  component  with  increasing  distance  from  the  plane. 
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Figure  3.20.  Wavenumber-frequency  plots  of  panel  vibration  at  kx  = 
0  (left)  and  ky  =  0  (right). 
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Figure  3.21.  Two-dimensional  wavenumber  spectrum  of  panel  vibra¬ 
tion  at  1000  Hz.  Note  that  the  scales  for  the  kx  and  ky  dimensions 
are  not  equal.  Note  that  the  x  and  y  axes  are  not  equal. 


The  half-monopole  source  is  therefore  thought  to  be  a  resonable  approximation 
of  a  monopole  for  testing.  While  the  evanescent  componenents  displayed  by  the 
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Figure  3.22.  Monopole  pressure  amplitude  at  the  distances  measured 
by  the  holography  array.  The  1/R  decay  in  pressure  is  clearly  evident. 


monopole  cannot  be  represent  a  broad  spectrum  of  complex  angle  components  rather 
than  a  single  plane  wave,  they  exhibit  similar  decay  over  distance,  properties  that 
allow  them  to  be  an  experimental  stand-in  for  evanescent  plane  waves. 

The  resulting  panel  vibration  due  to  acoustic  excitation  from  the  monopole  source 
is  shown  in  Figure  3.25  for  standoff  distances  of  22  and  202  millimeters.  The  22 
millimeter  standoff  shows  significant  vibration  in  wavenumber  components  above  the 
material  wavenumber  in  air.  While  supersonic  components  in  air  would  be  able 
to  excite  the  panel  vibration  due  to  the  presence  of  the  slow  flexural  waves,  such 
excitation  would  be  seen  across  a  range  of  standoff  distances,  as  the  energy  in  the 
supersonic  domain  does  not  decay  as  dramatically.  However,  the  202  millimeter 
standoff  does  not  exhibit  these  same  characteristics.  This  is  indicative  of  evanescent 
wave  propagation  into  the  material. 
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Figure  3.23.  Wavenumber  spectrum  of  the  monopole  sound  radia¬ 
tion  at  63  millimeter  offset.  The  assymetric  nature  of  the  wavenumber 
spectrum  is  due  to  the  offset  measurement  of  sound  pressure  radia¬ 
tion.  The  black  circle  indicates  the  material  wavenumber  in  air  of  18.3 
rad/m;  the  wavenumber  spectrum  is  clearly  highest  in  this  wavenum¬ 
ber  bin. 


3.5  Conclusions 

The  use  of  evanescent  waves  has  been  shown,  in  theory,  to  increase  energy  trans¬ 
mission  across  the  air-solid  interface.  By  representing  evanescent  waves  using  complex 
angles,  these  waves  can  be  input  into  classical  formulations  for  energy  transmission 
across  the  air-fluid  and  air-solid  interfaces.  The  evanescent  wave  will  overcome  the 
critical  angle  that  typically  limits  energy  propagation,  and  allow  for  transmission 
across  a  wider  range  of  angles. 

Beyond  the  critical  angle,  a  plane  wave’s  intensity  does  not  propagate  into  the  sys¬ 
tem,  existing  only  as  a  surface  pressure  that  rapidly  decays.  The  use  of  the  evanescent 
angle  does  not  increase  the  pressure  transmission  coefficient,  but  reduces  the  refracted 
angle,  allowing  for  propagation  of  energy  into  the  system.  The  use  of  complex  an¬ 
gles  introduces  an  imaginary  component  to  the  refracted  angle  of  the  propagating 
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Figure  3.24.  Measured  ratio  of  evanescent  to  non-evanescent  compo¬ 
nents  for  the  monopole  source  at  1000  Hertz. 
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Figure  3.25.  Wavenumber  spectrum  of  monopole  excitation  for  stand¬ 
off  distances  of  22  millimeters  (left)  and  202  millimeters  (right).  The 
decreasing  energy  in  the  system,  particularly  at  subsonic  wavenum¬ 
bers  (which  is  demarcated  by  the  dashed  black  line),  shows  the  effec¬ 
tiveness  of 
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wave;  this  will  cause  additional  decay  into  the  system  which  should  be  accounted  for 
by  measuring  the  intensity  propagating  into  the  system.  An  approach  for  intensity 
propagation  taking  into  account  the  pressure  transmission  and  the  resulting  intensity 
due  to  the  angle  of  transmission  must  be  formulated  to  allow  for  a  means  to  see  the 
effect  of  complex  angles  on  energy  transmission  across  the  interface. 

A  monopole’s  trace  pressure  distribution  along  a  surface  has  significant  evanes¬ 
cent  wave  components;  these  components  decrease  with  standoff  distance  between 
the  monopole  and  the  plane  (relative  to  the  supersonic  components  typically  seen). 
However,  a  monopole  source  used  to  excite  a  high-impedance  solid  object  will  exhibit 
these  evanescent  components  in  the  resulting  vibration  of  the  solid  object. 
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CHAPTER  4.  GENERATION  OF  NEAR-FIELD  PRESSURE 

DISTRIBUTIONS 

The  use  of  plane  waves  is  convenient  in  theoretical  formulations  of  fluid-solid  trans¬ 
mission  because  their  trace  pressure  along  a  plane  surface  is  a  steady-state,  fixed- 
wavelength  acoustical  interaction.  An  evanescent  wave,  formed  by  incorporating  an 
imaginary  component  into  the  angle  of  incidence,  will  have  a  fixed  decay  that  can  be 
used  in  theoretical  formulations. 

However,  plane  waves  are  impossible  to  generate  from  a  simple  sound  source  in  a 
free-held  environment.  While  plane  waves  can  be  generated  in  confined  spaces  such 
as  tubes,  or  in  certain  conditions  such  as  locally  in  the  far-held  of  a  monopole,  such 
cases  are  limited  by  the  geometries  and  frequencies  involved  in  the  system.  In  order 
to  utilize  typical  sound  sources  in  general  cases,  we  need  to  be  able  to  generate  plane 
waves  to  use  as  inputs  to  the  material  interface. 

The  focus  of  the  work  described  in  this  chapter  was  on  generating  plane  waves, 
particularly  evanescent  waves,  by  using  simple  sources  in  a  phased  array.  By  adapt¬ 
ing  a  previous  model  for  the  least-square  ht  of  pressure  distributions,  an  approximate 
pressure  distribution  can  be  created  by  the  phased  array.  An  exploration  of  the  factors 
affecting  the  generation  of  pressure  distributions,  including  the  incident  wave  param¬ 
eters,  source  geometry,  and  amount  of  decay  in  evanescent  waves,  will  be  explored 
here. 

The  model  developed  and  analyzed  here  is  for  a  two-dimensional  sound  pressure 
distribution  and  array.  While  this  does  not  provide  a  simulation  of  a  true  physical 
setup,  it  will  provide  an  overview  of  the  geometric  properties  of  the  array  and  char¬ 
acteristics  of  incident  waves  to  produce  pressure  distributions  on  the  surface.  For  a 
simulation  used  in  construction  of  an  experimental  array,  a  three-dimensional  model 
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should  be  created,  which  would  incorporate  the  polar  angle  of  the  incident  wave  as 
well  as  a  length  and  number  of  sources  along  the  additional  coordinate.  However, 
construction  of  an  experimental  array  should  also  take  into  account  the  possibility  of 
additional  geometries  and  the  angular  dependence  of  physical  speakers.  The  model 
here  is  two-dimensional  to  reduce  the  complexity  and  allow  for  ease  of  analysis  of  the 
simpler  parameters. 

4.1  Properties  of  Waves  and  Sources 

In  order  to  approximate  a  plane  wave  incident  upon  a  surface  and  maintain  its 
relationships  with  a  material  interface,  both  the  pressure  and  normal  velocity  of  the 
plane  wave  need  to  be  preserved  in  the  generated  pressure  distribution.  Accurately 
modeling  these  states,  as  well  as  the  pressure  generated  by  monopole  and  dipole 
sources,  for  example,  will  allow  for  the  development  of  a  model  to  approximate  plane 
wave  pressure  distributions  created  by  using  multiple  monopoles  or  dipoles. 

4.1.1  Properties  of  Plane  Waves 

The  properties  of  a  plane  wave  that  define  the  interaction  at  a  material  interface 
are  the  pressure  at  the  interface  and  the  velocity  normal  to  the  interface.  These 
quantities  must  be  constructed  on  the  surface  so  as  to  simulate  plane-wave  interaction 
with  the  surface. 

A  plane  wave  incident  on  a  surface  at  an  angle  9  has  a  pressure  distribution  in 
x  —  z  space  of: 


p  =  Aejk*x~jk*z  ,  (4.1) 

where  A  is  the  pressure  amplitude  and  kx  =  ksm(9 )  and  kz  =  kcos(9 )  are  the 
directional  wavenumbers  derived  from  the  material  wavenumber  k.  The  resulting 
normal  velocity  can  be  derived  from  the  pressure-velocity  relationship: 
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vz  =  —  ^z  =  Aejk*x~jk°zz  =  -C^lAejk*x~jkzZz  .  (4.2) 

jup  oz  up  cp 

For  the  purposes  of  our  model,  z  —  0  and  the  variation  of  p  and  vz  will  be  only 
dependent  on  x. 

If  the  trace  pressure  distribution  is  specihed  by  a  frequency  and  angle  of  incidence, 
the  resulting  normal  velocity  into  the  surface  is  unique  to  that  pressure  distribution. 
Therefore,  the  pressure  distribution  can  be  modeled,  and  the  velocity  should  match 
what  is  expected  from  an  incident  plane  wave  for  that  pressure  distribution. 

4.1.2  Properties  of  the  Acoustic  Monopole 

The  simplest  acoustical  source  is  the  monopole,  a  point  source  of  acoustic  pressure 
consisting  of  a  dilating  sphere  (of  diameter  much  smaller  than  a  wavelength)  with  a 
volume  output  of  Q  cubic  meters  per  second.  The  pressure  distribution  from  such  a 
monopole  is: 


kQ  e- 


,jkr 


^  c  (A  Q\ 

P=-JPoCo-A - ,  (4.3) 

47t  r 

where  p0  and  c0  are  the  density  and  speed  of  sound  in  air,  respectively,  k  is  the 
wavenumber,  and  r  is  the  distance  between  the  source  and  the  field  point.  The 
pressure  therefore  has  solely  a  radial  dependence,  and  decreases  from  the  origin  by 
1  jr.  The  velocity  resulting  from  a  monopole  is: 


1  dP  ^  kQ  1  ^ejkr  ^ 

Vr  =  —  —j—r  =  -3-7-  1  -  XT r) - r 

jup  dr  47 r  jkr  r 

where  vr  is  oriented  along  the  radial  axis  from  the  monopole.  If  the  monopolc  is 

located  above  the  plane,  the  resulting  normal  velocity  into  the  plane  will  be  vz{x)  = 

—  |fy  |  cos(0)T  ,  where  (j)  is  the  angle  between  the  normal  to  the  plane  and  the  line 

between  the  monopolc  location  and  the  point  of  incidence. 

The  trace  pressure  distribution  along  a  plane  surface  resulting  from  a  monopole 

varies  significantly  depending  on  the  standoff  distance  between  the  plane  surface 
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and  the  source.  In  the  far-held  where  the  1  jr  decay  of  the  sources  is  small,  the 
monopole’s  pressure  distribution  looks  like  a  non-evanescent  plane  wave  -  the  trace 
pressure  distribution  has  a  wavelength  approximately  equal  to  the  wavelength  of  a 
plane  wave  in  air,  and  there  is  very  little  decay.  However,  in  the  near-held,  the  1  jr 
decay  creates  an  evanescent  pressure  distribution,  though  not  one  that  conforms  to  an 
exponential  decay  or  that  has  a  constant  wavelength.  The  effective  wavelength  also 
changes  along  the  surface,  becoming  shorter  with  increasing  distance  from  the  source. 
As  a  result,  a  monopole  outputs  a  spectrum  of  wavelengths  and  amplitudes  which 
vary  according  to  the  position  of  the  point  of  incidence  relative  to  the  monopole. 

4.1.3  Properties  of  the  Acoustic  Dipole 

An  acoustic  dipole  consists  of  two  monopoles,  180  degrees  out-of-phase,  located  in 
close  proximity.  While  two  such  monopolcs  occupying  the  same  location  would  result 
in  complete  cancellation,  the  monopoles  forming  a  dipole  are  a  short  distance  apart, 
with  the  line  between  the  two  dipoles  forming  the  “dipole  axis”.  The  effect  of  this 
formulation  is  to  produce  angular  variance  of  the  output  pressure,  with  a  minima  of 
zero  radiation  normal  to  the  dipole  axis. 

The  dipole  pressure,  as  calculated  from  an  approximation  of  two  monopole  sources 
seperated  by  a  distance  that  is  small  compared  to  a  wavelength,  is: 

l'~(  )  o  pjkr 

P=-poC0 -7 - cos(C)  ,  (4.5) 

47T  r 

where  (  is  the  angle  relative  to  the  dipole  axis  and  s  is  the  distance  between  the 
two  dipoles  (each  monopole  in  the  dipole  has  a  source  with  amplitude  Q).  In  this 
formulation,  the  dipole  axis  will  be  oriented  normal  to  the  plane  containing  the  pres¬ 
sure  distribution.  For  the  purposes  of  these  calculations,  s  will  be  equal  to  1  /  10th 
of  a  wavelength  or  less.  The  output  of  a  dipole  is  cylindrically  symmetric  about  the 
dipole  axis,  with  positive  and  negative  pressure  regions  to  the  front  and  back.  The 
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maximum  radiated  pressure  is  along  the  dipole  axis,  reducing  with  angle  to  zero  at 
the  normal  to  the  dipole  axis. 

The  velocity  created  by  the  dipole  is  calculated  in  the  same  fashion  as  Equation 
4.4,  and  shows  the  same  relationship  to  the  monopole’s  velocity  as  the  monopole- 
dipole  pressure  relationship:  i.e., 

L‘~  ( )  Q  1  pjkr 

Vr  =  1  -  — )— cos(C)r  .  (4.6) 

47T  jkr  r 

Because  a  dipole  is  dependent  on  the  cancellation  between  two  monopole  sources, 
the  efficiency  of  the  dipole  radiation  is  much  lower  than  that  of  a  monopole  by  a  factor 
of  sk  cos(C);  for  s  equal  to  1  /10th  of  a  wavelength,  this  number  has  a  maximum  value 
of  0.63  and  the  dipole  will  always  be  a  less-efficient  radiator  than  a  monopole. 

The  dipole  is  useful  for  analysis  because  it  is  a  closer  approximation  of  a  real 
source  than  a  monopole.  An  unbaffled  speaker,  for  example,  will  exhibit  dipole-like 
behavior  clue  to  the  alternating  dilation  and  expansion  of  air  on  each  side  of  the 
speaker. 

While  a  single  monopole  or  dipole  cannot  generate  plane  waves,  a  phased  combina¬ 
tion  of  sources  can  be  used  to  approximate  an  incident  wave  on  a  surface.  By  adding 
pressure  from  multiple  sources,  we  can  approximate  a  desired  pressure  distribution. 

4.2  Least-Squares  Approximation  of  Pressure  Distributions  using  Multi¬ 
ple  Sources 

By  phasing  and  amplifying  each  source  in  an  array  appropriately,  the  array  can 
be  used  to  create  an  approximation  of  a  plane  wave  on  a  surface.  The  quality  of  the 
fit  between  the  ideal  and  generated  pressure  distributions  depends  on  many  variables, 
including  the  position  and  spacing  of  the  sources,  the  type  of  sources  used,  and  the 
amount  of  decay  in  an  evanescent  wave. 

A  schematic  of  the  situation  modeled  here  is  shown  in  Figure  4.1.  A  pressure 
distribution  along  a  plane  is  sampled  along  the  surface  over  a  distance  Lx  (1,024  points 
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were  sampled  in  this  model),  thus  creating  an  array  of  points,  P(xi),  representing  an 
incident  wave  of  a  particular  frequency,  /,  and  angle  of  incidence,  9r+j6i.  The  source 
array  is  centered  over  the  center  of  the  target  distribution,  with  N  sources  along  a 
line  of  length  L  located  at  Z0  above  the  plane. 


Ls 


N  (monopoles  or  dipoles) 


Figure  4.1.  Schematic  of  multi-source  array  used  to  fit  pressure  dis¬ 
tributions,  showing  variables  Lx ,  /,  6r  +  jOi,  Z0,  Ls,  and  N. 


The  method  for  generating  pressure  distributions  from  multiple  monopole  sources 
that  is  used  as  the  basis  for  this  calculation  was  first  presented  by  Kirkeby  and  Nelson 
[46].  Their  formulation  uses  monopole  sources;  expansion  of  this  model  to  allow  for 
dipoles  sources  may  provide  a  more  efficient  fit  to  the  desired  pressure  distribution. 
The  Kirkeby  and  Nelson  model  comprises  a  system  of  linear  equations  of  the  form: 

n  =  [HIM  .  (4-7) 

where  qj  is  the  complex  amplitude  (which  incorporates  a  phase  as  well  as  a  magnitude) 
of  the  source  located  at  Xj,  Pi  is  the  value  of  the  desired  pressure  distribution  at  points 
yi ,  and  [H]  is  a  matrix  that  relates  the  pressure  P(xi )  and  the  monopole  strength  qj. 
The  elements  of  the  matrix  are: 
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ejkrij 

Hij  =  - ;  r ij  =  | Xi  -yj  |  .  (4.8) 

rij 

The  solution  for  [q]  yields  the  monopole  phase  and  amplitude  for  a  least-squared- 
error  with  the  pressure  distribution.  The  resulting  pressure  approximation  could  be 
calculated  from  the  [q]  vector  by  summing  the  resulting  pressure  of  all  the  monopoles. 
The  monopole  source  strength.  Q ,  as  used  in  Equations  4.3  and  4.4  is: 


.  4vr 

°'=Wb 


(4.9) 


For  a  dipole,  a  similar  formulation  is  possible,  with  a  slight  modification  in  terms 
of  [H]  to  reflect  the  radiation  properties  of  a  dipole:  i.e., 


—  — - cos(Cij)  ,  (4.10) 

rij 

where  Q3  is  the  angle  between  ay  and  y3  (i.e.,  the  polar  angle  from  the  dipole  axis). 
Once  qj  has  been  defined  as  the  monopole  amplitude,  the  source  strength,  Qj,  of  the 
dipole  is: 


n  -  ■  4n 
J  J pocosk2^  ' 

In  order  to  determine  the  degree  to  which  the  generated  pressure  and  velocity 
distributions  fit  the  desired  incident  wave,  fit  coefficients  Rp  and  Ry  are  calculated 
from  the  sum  of  the  residual  values  at  each  ay  point  between  the  desired  pressure  P 
(or  V  for  Ry)  and  the  approximate  pressure  P  and  the  RMS  value  of  the  signal: 


Rp  =  1  — 


Ei(P(xi)  -  P(xi))2 


(4.12) 


E*  p(xi)2 

This  fit  coefficient  value  is  different  from  the  coefficient  of  determination  in  that  the 
sum  of  the  squared  value  is  always  taken  with  respect  to  zero,  the  sum  of  a  complete 
period  of  a  sine  wave.  For  an  oscillating  value,  the  sum  of  the  squares  will  yield  a 
value  of  Rp  that  has  a  maximum  value  of  1,  indicating  that  P(ay )  =  P(xi)  at  all  ay 
values.  A  value  of  0.98  or  higher  would  indicate  a  good  fit  of  the  pressure  distribution. 
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A  value  of  0  would  be  the  theoretical  minimum  value  obtained  during  a  least-squares 
curve  fit. 

Figure  4.2  shows  example  pressure  distribution  fits  for  a  generic  trace  pressure 
distribution.  While  the  fit  coefficient  0.883  approximates  the  pressure  distribution 
well,  the  0.999  fit  is  indistinguishable  from  the  ideal  pressure;  the  maximum  difference 
among  the  1024  point  in  that  case  is  0.0015,  less  than  0.2%  of  the  wave  amplitude. 


Figure  4.2.  Sample  approximate  pressure  distribution  for  a  plane 
wave  pressure  distribution;  approximate  pressure  distributions  with  a 
fit  of  0.883  and  0.999  are  shown. 

In  simulating  incident  plane  waves  and  the  required  multi-source  generation,  seven 
variables  come  into  play.  The  length  of  the  pressure  distribution,  Lx,  will  affect  the 
distance  that  needs  to  be  modeled,  as  well  as  accounting  for  the  variance  of  in  the 
overall  signal.  The  frequency,  /,  and  angle,  0r+j9i ,  of  the  incident  wave  will  determine 
the  surface  wavenumber,  with  higher  wavenumbers  having  more  variation  along  the 
surface  and  therefore  being  harder  to  approximate.  An  imaginary  component  of  the 
incident  angle  also  determines  the  relationship  between  pressure  and  velocity.  There 
are  N  sources  located  at  a  height  Z0  above  the  plane,  and  along  a  line  of  length  Ls 
centered  above  the  pressure  distribution.  Table  4.1  summarizes  these  variables.  Each 


59 


of  these  variables  affects  the  resulting  fit  of  the  equation  (expressed  by  Rp  and  Ry), 
and  should  be  considered. 

Table  4.1.  Table  of  variables  applicable  to  multi-source  arrays  and 
the  desired  pressure  distribution. 


Variable  Name 

Variable 

LInits 

Length  of  Distribution 

Lx 

meters 

Frequency  of  Wave 

f 

Hertz 

Angle  of  Incidence  (Real) 

dy 

radians 

Angle  of  Incidence  (Imaginary) 

0i 

radians 

Number  of  Sources 

N 

(none) 

Array  Standoff 

Z o 

meters 

Array  Length 

Ls 

meters 

4.3  Input  Parameter  Effects  on  Plane  Wave  Fit  Using  Multi-Source  Ar¬ 
rays 

Variations  in  the  seven  variables  determining  the  fit  of  an  incident  plane  wave 
using  a  multiple  source  array  effect  the  fit  quality  of  the  resulting  pressure  and  velocity 
approximations.  By  restricting  these  variables  to  the  important  factors  and  relevant 
values  it  is  more  easily  possible  to  analyze  a  complex  problem. 

In  order  to  reduce  the  number  of  variables  modeled,  we  will  assume  that  the  pres¬ 
sure  distribution  to  be  modeled,  Lx,  is  always  50  centimeters  in  length.  Clearly, 
greater  lengths  will  be  more  difficult  to  fit  (depending  on  the  complexity  of  the 
pressure  distribution),  requiring  more  sources  or  a  longer  array;  however,  this  is  a 
reasonable  approximation  for  the  size  of  a  suspicious  package. 

Since  the  interest  in  our  study  is  the  linear  region  of  wave  propagation  in  air,  our 
analysis  will  be  limited  to  frequencies  below  20  KHz.  The  incident  angle  can  vary 
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between  0  and  90  degrees;  however,  previous  analysis  has  shown  that  the  significant 
contributions  of  evanescent  waves  lie  at  lower  angles  of  incidence. 

The  number  of  sources  used  will  vary  from  1  to  64;  64  channels  would  be  allow 
8x8-channel  output  cards  to  generate  the  array  signals.  It  is  unknown  what  stand¬ 
off  distance  will  be  requested;  a  maximum  of  100  meters  is  assumed.  Similarly,  a 
maximum  source  length  of  15  meters  will  be  used. 

4.3.1  Incident  Wave  Parameters 

Of  primary  importance  is  whether  or  not  the  trace  pressure  distribution  generated 
is  properly  modeling  the  normal  velocity  inherent  in  the  incident  plane  wave.  The 
trace  pressure  distribution  is  a  unique  function  of  the  angle  of  incidence  and  frequency, 
and  the  appropriate  velocity  should  follow  if  the  pressure  is  fit  correctly.  However, 
approximating  the  trace  pressure  from  multiple  sources  will  not  necessarily  yield  the 
appropriate  velocity  relationship. 

For  a  plane  wave  distribution  generated  by  64  sources  with  a  source  offset  5  meters 
from  the  plane  and  a  span  of  4  meters,  the  resulting  fit  coefficient  variation  with  angle 
and  frequency  is  given  in  Figure  4.3.  While  the  pressure  has  a  gradual  transition  in 
the  quality  of  the  fit,  the  velocity  transitions  suddenly,  with  the  good  fit  region  of 
velocity  matching  the  good  fit  of  the  pressure.  This  correlation  between  matching 
pressure  and  velocity  at  better  regions  of  fit  shows  that  fitting  the  pressure  will  allow 
for  accurate  approximation  of  the  normal  velocity.  However,  poor  fit  of  the  pressure 
will  led  to  a  lower-quality  fit  of  the  velocity.  Therefore,  only  a  high  quality  of  fit  of 
the  desired  pressure  will  produce  the  proper  chracteristics  of  an  incident  wave. 

Figure  4.3  also  shows  the  effect  of  frequency  and  incidence  angle  on  the  ability  of 
the  multi-source  array  to  fit  a  desired  pressure  distribution.  An  increase  of  either  the 
frequency  or  angle  of  incidence  will  increase  the  variation  across  the  desired  region, 
thus  making  a  fit  more  difficult.  The  border  constituting  the  region  of  good  fit 
coefficient  appears  to  be  defined  by  an  inverse  relationship  between  the  sine  of  the 
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incident  angle  and  the  frequency,  sin(6l)  oc  1//.  This  would  seem  to  indicate  that 
the  trace  wavelength  Xx  =  c/(/sin(61))  could  serve  as  an  indicator  of  correct  fit; 
however,  Figure  4.4  shows  the  border  plotted  against  a  line  of  constant  wavelength 
fit.  The  wavelength  plotted  has  been  chosen  to  minimize  the  least-squares  error 
across  the  entire  frequency  range;  it  can  be  clearly  seen  that  the  border  of  good  fit 
does  not  follow  the  constant  wavelength  curve  precisely.  Because  the  wavelength  of 
the  monopole  varies  with  /,  the  fit  coefficient  is  not  solely  dependent  on  the  trace 
wavelength.  For  example,  a  1  kilohertz  wave  with  an  incident  angle  of  20  degrees 
(A*  =  1  meter)  will  not  have  similar  fit  coefficient  properties  (with  respect  to  the 
geometry  of  the  system)  as  a  wave  of  500  Hertz  with  an  incident  angle  of  43  degrees 
(which  has  the  same  equivalent  wavelength),  because  the  resulting  wavelength  of 
radiated  sound  from  the  monopole  sources  at  1  kilohertz  are  half  the  length  of  those 
radiated  by  500  Hertz  sources. 
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Figure  4.3.  Fit  coefficient  variation  for  the  pressure  (left)  and  normal 
velocity  (right)  with  angle  and  frequency  for  a  source  array  5  meters 
from  the  plane,  with  64  sources  spanning  4  meters. 


An  increase  of  the  length  of  the  source  array  will  increase  the  quality  of  the  fit, 
as  seen  in  Figure  4.5  for  the  geometry  of  Figure  4.3  except  with  a  source  span  of  10 


62 


Figure  4.4.  Border  of  R  —  0.98  fit  plotted  with  best-fit  line  of  trace 
wavelength,  \x  =  c/(/sin(0))  =  0.135.  The  fit  of  the  trace  wave¬ 
length  curve  shows  that  the  trace  wavelength  is  only  an  approximate 
indicator  of  the  good  fit  region. 


meters.  The  increase  in  the  quality  of  fit  is  possibly  due  to  the  additional  sources 
further  out,  which  will  produce  more  plane- wave-like  behavior.  The  additional  length 
of  the  source  array  reduces  the  source  density,  but  provides  improved  fit  at  higher 
angles,  which  have  a  lower  trace  wavelength.  This  would  indicate  that  there  isn’t  a 
proportional  relationship  between  source  density  and  the  trace  wavelength  that  can 
be  fit  by  the  array. 

Increasing  the  standoff  distance  will  reduce  the  quality  of  the  fit;  Figure  4.6  shows 
the  fit  for  a  source  array  with  a  length  of  10  meters  and  a  standoff  of  30  meters.  The 
same  relationsip  between  the  incident  angle  and  frequency  are  present;  however,  the 
fit  coefficients  are  much  lower  for  similar  angles  and  frequencies.  The  ’’spikes”  in  angle 
for  good  fit,  which  were  present  but  not  as  prominent  in  the  previous  simulations,  are 
due  to  the  poor  fit  of  non-integer  wavelengths  in  the  simulation.  An  increase  of  the 
pressure  distribution  length  can  shift  these  if  needed. 


Angle  (degrees)  Angle  (degrees) 
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Figure  4.5.  Fit  coefficient  variation  for  the  pressure  (left)  and  normal 
velocity  (right)  with  angle  and  frequency  for  a  source  array  5  meters 
from  the  plane,  with  64  sources  spanning  10  meters. 
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Figure  4.6.  Fit  coefficient  variation  for  the  pressure  (left)  and  normal 
velocity  (right)  with  angle  and  frequency  for  a  source  array  30  meters 
from  the  plane,  with  64  sources  spanning  10  meters. 


Fit  Coefficient  Fit  Coefficient 
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The  plots  of  fit  coefficient  show  the  difficulty  in  fitting  high  frequency  and  high 
angle-of-incidence  waves.  While  additional  geometry  and  source  parameters  may 
improve  the  fit  of  pressure  and  velocity,  the  angle  of  incidence  is  a  limiting  parameter. 

4.3.2  Geometry  Parameters 

As  seen  in  the  previous  section,  by  decreasing  the  standoff  distance  or  increasing 
the  source  length,  the  quality  of  the  fit  can  be  improved.  Fit  coefficient  variation 
with  these  parameters  is  shown  in  Figure  4.7  for  an  incident  wave  of  2  kilohertz  at  a 
20  degree  incidence  angle,  as  modeled  using  64  sources.  There  is  a  linear  relationship 
between  the  standoff  distance  and  the  source  length  for  the  border  of  fit  coefficients, 
although  significantly  more  source  length  is  required  to  obtain  an  improved  fit. 


5  10  15  5  10  15 
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Figure  4.7.  Fit  coefficient  variation  for  the  pressure  (left)  and  normal 
velocity  (right)  with  standoff  and  source  length  for  an  incident  wave 
of  2  kilohertz  at  an  incident  angle  of  20  degrees  modeled  using  64 
sources. 

For  a  higher  frequencies,  a  substantially  longer  source  length  must  be  used  to 
approximate  a  plane  wave,  as  seen  in  Figure  4.8.  The  lower  slope  of  the  good  fit 
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border  equates  the  to  need  for  a  greater  source  length  at  increased  standoffs.  For  a 
higher  incidence  angle,  the  quality  of  the  fit  decreases  further,  as  seen  in  Figure  4.9. 


5  10  15  5  10  15 
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Figure  4.8.  Fit  coefficient  variation  for  the  pressure  (left)  and  normal 
velocity  (right)  with  standoff  and  source  length  for  an  incident  wave 
of  5  kilohertz  at  an  incident  angle  of  20  degrees  modeled  using  64 
sources. 


A  study  of  the  source  geometry  shows  that  source  length  and  standoff  distance 
need  to  increase  in  tandem  in  order  to  maintain  a  good  fit  of  the  pressure  distri¬ 
bution.  For  higher  frequencies  and  incidence  angles,  the  proportional  source  length 
increase  with  standoff  becomes  greater;  that  is,  a  greater  source  length  is  needed  for 
an  equivalent  standoff  distance  in  harder-to-fit  plane  waves. 


4.3.3  Source  Length  and  Density 

Previous  examples  have  been  based  on  the  assumption  that  the  maximum  number 
of  sources  was  used;  a  reduction  of  the  number  of  sources  in  the  array  would  allow 
for  a  lower-power  solution  to  obtain  evanescent  waves. 
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Figure  4.9.  Fit  coefficient  variation  for  the  pressure  (left)  and  normal 
velocity  (right)  with  standoff  and  source  length  for  an  incident  wave 
of  5  kilohertz  at  an  incident  angle  of  40  degrees  modeled  using  64 
sources. 


Figure  4.10  shows  the  typical  variation  of  fit  coefficient  with  source  number  and 
length,  using  an  incident  wave  of  2  kilohertz  at  20  degrees  and  a  linear  array  30 
meters  from  the  surface.  The  region  of  good  fit  is  defined  by  a  horizontal  border 
on  the  lower  end  and  a  high-slope  vertical  boundary.  The  horizontal  slope  suggest 
a  minimum  number  of  sources  needed  to  approximate  the  incident  wave;  for  the 
current  situation,  this  is  8  sources,  far  fewer  than  the  64  in  modeling  efforts  with  a 
fixed  number  of  sources.  The  vertical  hue  shows  a  slight  slope,  indicating  a  decrease 
in  source  number  with  additional  length  required  to  maintain  fit. 

For  higher  frequencies,  as  shown  in  Figure  4.11,  the  minimum  number  of  sources  is 
increased;  this  would  suggest  that  the  number  of  sources  is  a  function  of  the  frequency 
and  incidence  angle.  There  is  a  small  slope  at  the  bottom  edge  of  the  good  fit  region, 
thus  showing  the  need  for  increased  source  numbers  with  additional  length;  this  is 
indicative  of  a  minimum  source  density  needed  to  accurately  model  the  incident  wave. 
Otherwise,  the  higher  frequency  plot  mirrors  the  trend  seen  in  the  previous  figure, 
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Figure  4.10.  Fit  coefficient  variation  for  the  pressure  (left)  and  nor¬ 
mal  velocity  (right)  with  source  length  and  number  of  sources  for  an 
incident  wave  of  2  kilohertz  at  an  incident  angle  of  20  degrees  and  a 
source  standoff  of  30  meters. 

where  this  is  very  little  variance  once  a  minimum  source  length  and  number  have 
been  achieved. 

4.3.4  Modeling  Evanescent  Waves 

By  treating  the  incident  angle  as  a  complex  angle  and  separating  its  real  and 
imaginary  components  into  separate  variables,  the  effect  of  evanescent  waves  can  be 
introduced  into  the  model.  The  use  of  complex  angles  will  introduce  a  pressure  decay 
into  the  system  and  a  phase  difference  between  the  pressure  and  velocity  that  are  not 
seen  when  modeling  conventional  plane  waves. 

The  addition  of  an  imaginary  component  to  the  angle,  which  is  shown  in  Figure 
4.12  for  variation  with  real  and  imaginary  angle  for  a  2  kilohertz  incident  plane  wave, 
has  only  a  slight  effect  on  the  fit  of  the  pressure,  reducing  the  angle  of  good  fit  by 
9  degrees  from  a  complex  angle  of  0  to  0.1  radians.  The  quality  of  the  velocity 
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Figure  4.11.  Fit  coefficient  variation  for  the  pressure  (left)  and  nor¬ 
mal  velocity  (right)  with  source  length  and  number  of  sources  for  an 
incident  wave  of  8  kilohertz  at  an  incident  angle  of  20  degrees  and  a 
source  standoff  of  10  meters. 


fit,  however,  decreases  significantly  with  increasing  evanescent  angle.  This  is  due 
to  the  phase  difference  that  the  evanescent  angle  imposes  on  the  velocity;  a  greater 
imaginary  component  pushes  the  phase  difference  between  the  pressure  and  velocity 
towards  90  degrees.  While  a  monopole  has  properties  similar  to  an  evanescent  wave 
in  its  near-held,  these  are  dependent  on  geometry,  and  relying  on  far-held  properties 
of  a  monopole  will  not  produce  evanescent  components. 

4.3.5  Dipole  Multi-Source  Arrays 

The  use  of  dipole  sources  in  the  array  has  little  effect  on  any  of  the  parameters  or 
the  quality  of  ht  shown  for  monopoles.  All  the  previously  presented  results  show  no 
appreciable  difference  when  modeled  using  dipoles  in  place  of  monopoles.  For  exam¬ 
ple,  Figure  4.13  shows  the  improvement  in  ht  coefficient  for  a  dipole  array  compared 
to  the  monopole  array  for  the  angle  and  frequency  variation  shown  in  Figure  4.5. 


69 


o.i 


0.09 


0,08 

^  0.07 
<o 

»  0.06 
u-. 

<  0.05 
> 

c  0.04 

a 

aJ 

E  0.03 


0.02 


0.01 


. 

20  40  60  80 

Real  Angle  (degrees) 


20  40  60  80 

Real  Angle  (degrees) 


Figure  4.12.  Fit  coefficient  variation  for  the  pressure  (left)  and  nor¬ 
mal  velocity  (right)  with  real  and  imaginary  angle  components  for  an 
incident  evanescent  wave  of  2  kilohertz  and  a  source  geometry  of  4 
meters  with  a  standoff  of  20  meters. 


While  there  are  small  regions  of  fit  with  an  increase  of  0.25,  all  these  regions  occur 
where  the  monopole  distributions  provide  a  poor  fit,  with  the  highest  fit  being  0.83. 

Figure  4.14  shows  the  improvement  in  fit  coefficient  when  using  a  dipole  array  for 
the  same  data  shown  in  Figure  4.10.  The  use  of  a  dipole  array  has  a  negligible  effect 
on  fit  coefficient,  with  a  maximum  improvement  of  0.01  in  small  regions,  and  change 
of  less  than  0.001  in  most  regions. 

That  the  dipole  and  monopole  arrays  have  such  close  fit  and  similar  parameters 
is  indicative  of  the  local  effect  of  each  individual  source.  A  dipole  differs  from  a 
monopole  by  the  cosine  of  the  angle  to  the  dipole  axis;  the  acoustical  output  near 
the  axis  will  be  similar  to  the  monopole.  Because  the  dipole  and  monopole  outputs 
from  the  array  are  similar,  this  suggests  that  the  effects  of  each  source  are  localized 
near  the  dipole  axis  (which,  in  this  case,  is  normal  to  the  plane).  There  may  be 
an  additional  benefit  to  using  dipole  sources  if  the  orientation  of  each  source  can  be 
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Figure  4.13.  Fit  coefficient  improved  between  the  dipole  and 
monopole  arrays  for  a  source  array  5  meters  from  the  plane,  with 
64  sources  spanning  10  meters. 
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Figure  4.14.  Fit  coefficient  improved  between  the  dipole  and 
monopole  arrays  for  an  incident  wave  of  2  kilohertz  at  an  incident 
angle  of  20  degrees  and  a  source  standoff  of  30  meters. 
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optimized  as  part  of  the  model;  however,  the  addition  of  dipole  orientation  to  the 
model  will  form  a  nonlinear  system  of  equations. 

A  dipole  is  a  much  less  efficient  radiator  than  a  monopole;  therefore,  for  an  equiv¬ 
alent  source  array  with  dipoles  that  have  the  same  pressure  output  as  monopoles,  the 
array  will  be  less  efficient  by  the  same  factor,  sk  cos(£),  that  a  dipole  is  less  efficient 
than  a  monopole. 

4.4  Conclusions 

Multi-source  arrays  have  been  shown  to  be  a  viable  option  for  constructing  trace 
pressure  distributions  on  a  plane.  The  model  constructed  here  shows  the  influence 
of  incident  wave,  source  geometry,  and  evanescent  wave  parameters  on  the  fit  of  the 
desired  waves,  for  both  monopoles  and  dipoles. 

While  monopoles  and  dipoles  have  characteristics  distinct  from  a  plane  wave,  a 
multi-source  array  can  be  modeled  using  a  least  squares  formulation  to  produce  an 
approximation  of  the  desired  incident  pressure  and  velocity.  The  quality  of  the  fit  of 
the  approximation  is  judged  by  comparing  the  desired  pressure  and  velocity  to  their 
approximations. 

The  fit  of  the  multi-source  system  is  heavily  dependent  on  the  incident  wave 
parameters;  beyond  an  angle  dependent  on  the  geometry  of  the  source  array,  the  fit 
is  poor,  with  further  and  smaller  arrays  being  able  to  fit  fewer  angles  of  incidence. 
Frequency  plays  a  role  as  well,  with  a  decrease  in  the  quality  of  the  fit  with  increasing 
frequency. 

In  the  geometry  of  the  source  array,  the  standoff  distance  decreases  the  quality  of 
the  fit,  but  additional  length  can  be  added  to  the  array  to  improve  the  quality  of  the 
fit.  There  is  a  minimum  source  density  that  must  be  achieved,  but  beyond  this  there 
is  little  variation  in  the  quality  of  the  fit. 

Evanescent  waves  adversely  affect  the  fit  of  the  system,  possibly  due  to  the  phase 
difference  between  the  pressure  and  velocity  in  the  system.  While  the  pressure  fit 
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varies  only  slightly  with  evanescent  components,  the  velocity  will  fit  will  be  worse 
due  to  phase  conditions  imposed  on  the  wave  that  cannot  easily  be  fulfilled. 

A  dipole  array  provides  similar  fit  as  a  monopole  array,  with  all  parameters  main- 
taing  similar  variation  to  their  monopole  counterparts. 
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CHAPTER  5.  WAVE  PROPAGATION  IN  MULTI-LAYER 

MATERIALS 

The  amplitude  and  phase  of  a  wave  propagating  through  an  infinite  homogeneous  ma¬ 
terial  change  only  due  to  the  phase  change  with  distance  and  possibly  the  damping 
in  the  material.  For  layered  materials  each  component  layer  has  its  own  propagation 
characteristics,  the  combination  of  which  causes  the  resulting  multi-layer  material  to 
have  a  distinct  susceptibility  to  acoustic  excitation.  By  modeling  the  wave  propaga¬ 
tion  characteristics  of  the  overall  material,  we  can  accurately  design  waveforms  that 
excite  the  components  of  the  material  that  are  of  interest. 

In  this  chapter,  a  model  for  wave  propagation  through  layered  materials  is  pre¬ 
sented  and  verified.  Applications  for  the  model  and  analysis  techniques  in  the  wavenumber- 
frequency  domain  are  presented.  A  discussion  of  the  effects  of  evanescent  waves  on 
layered  materials  concludes  the  chapter. 

5.1  Multi-Layer  Wave  Propagation  Model 

The  multi-layer  propagation  model  used  here  is  based  on  a  formulation  by  Brouard 
[21],  which  combines  the  previous  models  for  transfer  matrix  calculation  of  waves  in 
a  solid  (by  Folds  [19])  and  a  fluid  (Pierce  [20])  by  developing  these  transfer  matrices 
into  a  system  of  equations.  Another  major  innovation  in  Brouard’s  formulation  was 
the  use  of  a  fluid-solid  coupling  matrix,  which  would  transfer  the  material  states  in 
a  fluid  into  those  in  a  solid  (and  vice-versa).  In  these  models,  the  material  states  at 
each  interface  are  unknowns  and  propagation  characteristics  through  each  layer  are 
assumed  to  relate  states  on  each  side  of  a  layer  of  material.  The  transfer  matrices, 
which  Pierce  and  Folds  used  in  series  but  which  Brouard  uses  as  a  component  for 
developing  a  system  of  equations,  propagate  the  pressure  and  velocity  through  a  region 
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of  material  accounting  for  the  phase  change  over  a  length  L  in  the  material  layer.  For 
example,  the  relationship  between  the  pressure  and  normal  velocity  amplitude  vector 
[P]  in  a  fluid  from  a  Zq  to  zl  span  is  (for  an  e~^ut  time  convention): 


r  -| 

p 

cos(kzL)  jfsin(kzL) 

P 

(5.1) 

P 

= 

= 

z 

L  J 

z=0 

Vz 

z= 0 

j^sm(kzL)  cos  (kzL) 

Vz 

z=L 

where  kz  is  the  wavenumber  in  the  direction  normal  to  material  surface.  A  solid  ma¬ 
terial  will  have  an  equivalent  [P]  vector  containing  pressure  and  velocity  amplitudes 
for  both  longitudinal  (az,  vz)  and  shear  waves  (ax,  vx ),  and  a  4-by-4  transfer  matrix 
relating  pressure  and  velocity  relationships  over  distance. 

Transfer-matrix  methods  of  simulating  wave  propagation  through  multi-layer  ma¬ 
terials  relate  the  material  states  between  layer  1  and  layer  N  as  a  series  of  these 
matrices  that  reflect  propagation  through  each  layer: 


(5.2) 

The  Brouard  model  expands  upon  this  transfer- matrix  approach  to  develop  a  sys¬ 
tem  of  linear  equations  (a  similar  formulation  will  be  expanded  upon  in  the  following 
section): 
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The  resulting  matrix  is  an  N  —  1  by  N  system;  the  transfer  function  between  the 
incident  pressure  amplitude  p\  and  the  m-th  pressure  or  velocity  amplitude  in  the 
matrix  formulation  is  related  to  the  derivative  of  a  subset  of  the  matrix  with  the  first 
and  m-th  row  removed: 
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Hlm  =  (l  +  R1)^.  (5.4) 

l-^i  I 

While  this  type  of  formulation  allows  for  calculation  of  multiple  states  in  the  sys¬ 
tem,  the  drawbacks  are  that  it  does  not  distinguish  between  forward-  and  backward¬ 
going  waves  in  the  system,  and  that  inter-layer  material  states  cannot  be  discerned. 
Both  of  the  latter  quantities  are  important  for  calculating  the  intensity  transmitted 
through  layers,  which  would  be  a  true  indicator  of  energy  propagation  into  the  system. 
Furthermore,  the  use  of  a  wave-based  method  would  allow  for  evanescent  wave  inputs 
in  order  to  model  their  effectiveness  in  energy  propagation.  The  model  described  here 
modifies  a  formulation  by  Jessop  et  al.  [26],  which  uses  the  wave  potentials  in  each 
layer  to  simulate  propagation.  The  resulting  system  of  equations  for  the  multi-layer 
system  are  generated  by  equating  the  states  at  each  interface. 

5.1.1  Wave  Potential  Model  Formulation 

The  multi-layer  material  modeled  consists  of  (. N—2 )  layers  of  fluid  or  solid  bounded 
by  semi-infinite  fluid  layers  on  each  side,  with  the  incident  fluid  numbered  layer  1 
and  the  transmitted  bounding  fluid  layer  N.  An  incident  wave  in  the  semi-infinite 
first  layer  is  incident  at  an  angle  0\.  The  trace  wavenumber  on  the  surface  kx  = 
kix  —  A;1sin(d1)  is  constant  throughout  all  layers,  while  the  normal  wavenumber 
kiz  =  k±  cos (0i )  will  change  in  each  layer,  causing  refraction  in  each  subsequent  9t.  A 
schematic  of  the  setup  is  shown  in  Figure  5.1. 

The  system  of  equations  governing  wave  propagation  in  the  multi-layer  material  is 
formulated  by  equating  the  states  at  each  of  the  (IV  —  1)  interfaces  in  the  system.  The 
forward-  and  backward-going  pressure  and  shear  wave  amplitudes  are  four  unknowns 
in  each  layer;  for  a  fluid,  the  shear  wave  amplitudes  are  set  to  zero.  The  material 
states  are  calculated  as  follows. 

A  wave  propagating  in  the  negative  ^-direction  through  a  fluid  or  solid  at  an  angle 
6  can  be  represented  by  a  wave  potential  of: 
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Figure  5.1.  Schematic  diagram  of  a  sample  N  =  5  multi-layer  mate¬ 
rial.  A  semi-infinite  fluid  exists  on  each  side  of  the  three-component 
material. 


</>  =  Aejkxx-jkzz  .  (5.5) 

(The  time  constant  e~^ut  is  excluded  for  brevity).  Note  that  this  potential  is  similar 
to  the  pressure  equation  for  a  propagating  wave;  however,  the  amplitude  A  of  the 
wave  potential  is  not  equivalent  to  the  pressure  amplitude. 

A  single  wave  potential  equation  is  used  to  represent  longitudinal  (pressure)  waves, 
and  another  wave  potential  ^  (which  is  similar  to  the  equation  for  <fi  but  propagates 
at  an  angle  7  with  wavenumber  k)  is  used  to  represent  transverse  (shear)  waves  in 
the  system: 


(/,  _  jn~~- 


(5.6) 
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Figure  5.2.  Wave  propagation  characteristics  and  geometry,  show¬ 
ing  the  four  waves  in  a  solid  system:  two  forward-going  waves,  both 
longitudinal  and  shear,  and  their  companion  backward-going  waves. 


At  a  material  interface,  a  wave  incident  in  one  material  will  generate  a  reflected 
wave  in  the  incident  material  and  transmitted  waves  in  the  other  material.  The 
reflected  wave  travels  in  the  positive  z  direction  and  is  expressed  as: 

cj)  =  ARejk*x+jk°z  ,  (5.7) 

with  an  equivalent  expression  for  shear  waves  in  a  solid  material.  The  transmitted 
waves  travel  in  the  negative  z  direction  and  have  the  expected  form  with  amplitude 
Aj'. 

A  solid  material  will  be  defined  by  a  density  p  and  longitudinal  and  shear  wave 
speeds  c  and  b.  The  wave  speeds  can  be  calculated  from  other  material  properties; 
for  example,  if  given  the  Lame  constants  A  and  p,  the  wave  speeds  c  and  b  can  be 
calculated  as: 


'  A  +  2yU 


c  = 


P 


(5.8) 
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b=x[^.  (5.9) 

V  P 

The  stresses  in  a  solid  consist  of  the  normal  stress  az  and  tangential  stress  ax,  which 
can  be  computed  from  the  wave  potentials  as: 


<tz  =  ~ 

u 


.  .  ,  <920  <920  \ 


(5.10) 


_  .7 

&X 

U 


d2(j)  <920  <920 

I1  ^  dxdz  dx 2  dz2 


(5.11) 


A  fluid  cannot  sustain  shear;  therefore  0  —  0,  p  =  0,  and  ax  =  0,  and  the  pressure  in 
the  fluid  simplifies  to: 


p  =  —  AV20  =  —  icupcp  .  (5-12) 

The  velocity  in  fluid  or  solid  materials  can  be  computed  from  the  wave  potentials  as: 


v  =  V</>  +  V  x  -0  , 

which  can  be  simplified  into  its  constituent  parts: 


dd>  dib 

vz  =  w  b  vp  , 

az  ox 


(5.13) 


Ur  = 


<90  <90 


(5.14) 


(5.15) 


dx  dz 

In  a  multi-layer  material  the  model  is  formulated  by  using  all  the  boundary  conditions 
to  create  a  system  of  linear  equations  in  the  wave  amplitude  coefficients. 

Each  constrained  layer  (layers  2  through  N  —  1)  has  forward-going  (into  the  ma¬ 
terial,  traveling  in  the  negative  z  direction)  and  backward-going  (traveling  in  the 
positive  z  direction)  waves  that  are  the  product  of  transmission  and  reflection  at 
each  interface,  labeled  as  0 i+  and  0j_,  respectively  (with  equivalent  0,+  and  for 
shear  waves  in  a  fluid).  The  total  longitudinal  potential  0*  in  a  layer  represents  the 
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combined  effects  of  both  forward-  and  backward-going  waves;  the  material  states  in 
Equations  5.10-5.15  are  calculated  using  the  sum  of  both  waves.  The  wave  amplitudes 
are  represented  in  a  vector  P%  for  each  layer. 

At  the  interface  between  the  i-th  and  (i  +  l)-th  layers,  located  at  a  depth  of  the 
relationships  between  the  material  states  form  a  system  of  equations  that  constitute 
several  rows  of  the  system  of  equations  for  the  material.  For  solid-solid  or  fluid- 
fluid  layer  interfaces,  the  pressure  and  velocity  in  each  direction  is  equal  across  each 
interface,  and  the  pressure  and  velocity  amplitude  [P]i  at  the  right  side  of  the  i-th 
layer  is  equal  to  the  pressure  and  velocity  at  the  left  side  of  the  ( %  +  l)-th  layer. 

For  example,  the  states  across  a  fluid-fluid  interface  can  be  equated  by  applying 
equations  5.10  and  5.14  to  yield  two  equations: 


jupi4>i(zi )  -  (- jupi+1)(j)i+1(zi )  =  0  , 


(5.16) 


5(pi(zi)  8(pi+1(zi ) 


=  0 


(5.17) 


8z  5z 

Substituting  the  wave  potential  equations  into  5.16  for  the  forward-  and  backward¬ 
going  waves  gives: 


-juipi(Al+e  lk‘-m  +  =  -juipl+1(Am)+e  +  A{l+I)_e1Kj*m)  , 

(5.18) 

and  equation  5.17  will  simplify  to: 


+  {jkZti)Ai_ejkz'iZi  = 
-jkz,i+1A{i+1)+e-jk*’i+lZi  +  jkZii+lA{i+l)_ejk*’i+lZi  =  0  . 


(5.19) 


The  resulting  system  of  equations  for  the  interface  are: 
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Ai+ 

— pie  jkz’iZi 

—pieikz’iZi  Pi+ie  Az,i+izi 

pi+iejkz’i+lZi 

Ai- 

0 

—kz^e~^kz’iZi 

kZjiejkz’iZi  kZ:i+1e~jkz’i+lZi 

-kZji+1ejkz’i+lZi 

^(1+1)+ 

0 

^q+i)- 

(5.20) 

For  a  solid-solid  interface,  the  four  equations  for  continuous  longitudinal  and  shear 
stress  and  velocity  form  four  equations: 


(M-^i  -  *£i)  -  2 ^kl)(Al+e-^  +  A.e^) 
+2 HiKXiiKZii(Bi+e~iKz'iZi  -  Bi_ejKz’iZi)  = 
(\i+i(-kli+1  -  ei+1)  -  2 iii+1kli+1)(Aii+1)+e~jk^Zi  +  A(i+1)_ejkz^Zi) 
+2ni+1KXti+1KZji+1{B(i+1)+e-jK^+lZi  -  B(i+1)_ejK*’i+lZi) 
On. Jc  Jc  -(A-  p—foz,iZi  _  A.  pjkz,i%i\ 

+^(-<i  +  ^z,i){Bi+e-j^iZi  +  Bi_ejK*’iZi)  = 

2fii+1kXti+1kZ:i+1(A{i+1)+e-jkz-i+lZ'  -  A{i+1)_ejk*’i+lZi ) 
+/ii+i(-<i+1  +  <i+1)(£(i+ l)+e-jK^+lZi  +  fi(i+1)_e^z'i+1^) 
jkz,(-Ai+e~jk^  +  A_ejk^) 

+jKx,i(Bi+e~jKz’iZi  +  Bi_ejKz’iZi)  = 
jkZii+l(-A{i+1)+e-jk°’i+lZi  +  A(i+i)_eJ'fc*-i+1*i) 

+J«a;,i+i(5q+i)+e-i^'i+1^  +  B(i+i)_eiKz'i+1^) 
u  Y4.  p-jkz,iZi  ,  /i  pjkZ}izi\ 

+jKz,i(Bi+e-jK*’iZi  -  11,  f  )  = 

+J^,*+i(5(i+i)+e-iKz'i+121  -  fi(i+i)_ejKz’i+1^)  , 


(5.21a) 


(5.21b) 


(5.21c) 


(5.21d) 


which  can  be  easily  manipulated  into  a  system  of  linear  equations  similar  to  5.20. 
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At  a  fluid-solid  interface,  the  normal  stresses  (equivalent  to  pressure  in  the  fluid) 
and  normal  velocities  are  equal.  The  fluid  layer  cannot  sustain  shear  stress,  so  the 
shear  stress  in  the  solid  at  the  interface  is  zero.  There  is  no  equation  relating  the 
shear  velocity  at  the  interface.  The  resulting  material  states  at  an  i-th  layer  of  fluid 
and  an  (i  +  l)-th  layer  of  solid  are: 

Pi  &z,i+ 1  j  0  @x,i+ 1  j  Vi  hz,i+l  •  (5.22) 

and  the  resulting  equations  will  be: 


jupi(Ai+e  jkz-iZi  +  Ai_ejkz'iZi )  = 

-((Xi+i(-k2x,i+i  ~  k%+ 1)  -  2ni+lkli+l)(A{i+1]+e-jkz^lZi  +  A(i+1)_ejkz^Zi)  (5.23a) 

id 

+2pi+1(KXti+1KZti+1)(Bii+1)+e-jKz’i+lZi  -  B(i+1)_ejK^+lZi)) 

0  =  2pi+12kx,i+1kz,i+1(A{i+l)+e~jkz’i+lZi  -  A{i+1)_ejkz’i+lZi) 

(5.23b) 

+/h+i(-<*+i  +  1)+e~jK^+lZi  +  B{i+1)_ejKz^Zi) 

jkZji(—Ai+e~jkz’iZi  +  A_e^Zi)  = 

jkz,i+1(-A{t+1)+e-jkz^z'  +  A(i+1)_ejkz^Zi)  (5.23c) 

+jKx,i+1(B(i+1)+e-jKz’i+lZi  +  B(i+1)_ejKz’'+1A  ■ 


The  interface  between  the  first  two  layers  involves  the  incoming  wave  with  unit 
amplitude;  the  resulting  equations  differ  from  other  interface  equations  in  that  they 
form  an  inhomogeneous  system  of  equations  due  to  the  forcing  of  the  incident  wave. 
Only  the  reflected  wave  is  an  unknown  coefficient  in  the  first  layer;  the  amplitude  of 
the  incident  wave  is  a  constant.  For  example,  at  a  fluid-fluid  interface  between  the 
incident  layer  and  the  second  layer  at  z  —  0,  the  resulting  state  equations  will  be: 


-  jupiAi_  -  {-jupi{A2+  +  A2-))  =  -(-jupi)  ,  (5-24) 


jkZtiAi_  -  (- jkZt2A2+  +jkZt2A2_)  =  —(—jkZji)  . 


(5.25) 
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On  the  final  layer  N  (which  must  be  a  fluid  in  order  to  create  a  consistent  system 
of  linear  equations),  only  the  outgoing  wave  is  allowed  for.  Because  the  fluid  is  semi¬ 
infinite,  there  will  be  no  backward-traveling  component. 

The  systems  of  equations  for  each  interface  can  be  summarized  as  the  relationships 
between  coefficients  for  each  layer  [S'*],  a  vector  of  wave  amplitudes  [P,] ,  and  the 
resulting  constants  [C(]  (which  will  be  zero  for  all  interfaces  except  the  first  one): 

Si]  [p]  -  [si+i]  [pi+ 1]  =  -  [Ci]  ,  (5.26) 

and  these  constituent  parts  for  each  layer  can  be  combined  into  a  multi-layer  system 
of  the  form: 


(5.27) 

(5.28) 

or,  in  simplified  form: 

S  P  =  C  ,  (5.29) 

and  the  resulting  coefficients  [P]  can  be  solved  as: 

P  =  S  C  ■  (5.30) 


Once  the  coefficients  in  the  system  have  been  found,  the  amplitudes  of  each  wave 
can  be  input  into  their  respective  equations,  and  the  steady-state  material  states  can 
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be  found  at  all  points  z  by  summing  the  results  from  the  forward-  and  backward- 
going  waves  in  the  appropriate  layer.  These  material  states  are  transfer  functions 
between  the  state  and  the  unit  input  wave  amplitude;  to  obtain  a  transfer  function 
for  incoming  pressure,  the  material  states  should  be  divided  by  —;jujp\. 

Intensity  in  the  system  can  be  found  from  the  stress  and  velocity  (we  will  use  the 
generic  cr  and  v  to  correspond  to  either  the  longitudinal  or  shear  values)  at  any  state: 

J(-)  =  ,  (5.31) 

where  v*(z)  is  the  complex  conjugate  of  the  velocity  v(z).  In  a  fluid,  the  inten¬ 
sity  in  the  i-th  layer  (the  i  subscript  will  be  excluded  from  material  properties  and 
wavenumbers  for  clarity;  note  that  z  is  a  position)  is: 

I{z)  =  ^R(pu*z)  =  ]^t(-jup(A+e-jk*z  +  A_ejk*z)(jkz{-A+e-jk*z  +  A_ejk*z)*)  = 
^3?(- jup(A+e~jkzZ  +  A_ejk*z)(-jk*z(-A*+ejk*z  +  A*_e~jkzZ))  = 
K(-A+A*+  -  A_A*+e2jk*z  +  A+A*_e~2jkzZ  +  A_A*_ ))  . 

(5.32) 

The  center  two  terms,  which  contain  phase  components,  create  a  purely  imaginary 
term  that  will  not  contribute  to  the  intensity;  therefore,  intensity  is  constant  across 
a  fluid  layer. 

In  a  solid  layer,  the  intensities  in  the  longitudinal  and  shear  waves  both  consists 
of  16  terms,  obtained  from  multiplying  the  4  terms  constituting  stress  by  the  four 
terms  constituting  velocity.  The  longitudinal  intensity  is: 


84 


I,(z)  =  ht(i 

Z  oJ 

((A {-k2x  -  k\ )  -  2fik2z)A+e~jk*z  +  2 \fjiKxKzB+e-jK’z+ 
(A (~k2x  -  k2z)  -  2 /ik2z)A_ejk*z  -  2iiKxKzB_ejK*z) 
(jk*zA*+ejkzZ  -  jn*xB*+ejKzZ  -  jk*zA*_e~jkzZ  -  jK*xB*_e~jKzZ))  = 

2-»(-K(\(-kl  -  kl)  -  2^)A+Al 
+<(A (-kl  -  k\)  -  2„kl)A+B+’e*-*‘+^‘ 
+K(K~kl  -  kl)  -  2 pkl)A+A-  e-*k-‘ 

+<(K-K  -  kl)  -  2 

~k*z2imxKzB+A+*  eJ'(-«*+fc*)* 
-\-k*x2iikxkzB+B+ 
+k*z2nKxKzB+A_*  ej(-K°-k°> 
+k*x2  hkxkzB+B_*  e~2jKzZ 
-K{\ {-kl  -  kl)  -  2 pk2z)A_A+*e2^z 
+<(\(-kl  -  kl)  -  2vkl)A_B+’e<V“+-> 
+K(M-kl  -  kl)  -  2nkl)A_A_' 
+<(A (-fcf  -  k2z)  - 

+kt2iiK,KzB-A+ e!(*’+k’)* 
—k*x2iikxkzB_B+  e2^KzZ 
-k*z2nKxKzB_A_*  ej(Kz~kz> 
—  K*x2llKxKzB-B-  )  , 


(5.33) 
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while  the  shear  intensity  is: 

4  (z)  = 

Z  uj 

(2  nkxkzA+e~jk*z  +  h{-kI  +  K2z)B+e-jK’* 
-2  nkxkzA_ejk*z  +  h(-k2x  +  K2g)B-e?*'z) 
{-jk*xA*+ejk°z  -  jK*zB*+ejK*z 
-jk*xA*_e~jkzZ  +  jK*zB*_e~jK*z))  = 

2  fik*xkxkzA+A*+ 
-2  K*znkxkzA+B*+ej(~k*+K*> 
+k*x2nkxkzA+A*_e~2jk*z 
-2  K*zixkxkzA+B*_e~j{-kz+Kz)z 

+KK-*l  + 

+**//(-«£  +  ^)5+hHe-i(Kz+fcz)z 
+  4)B+B*_e~2^z 
-k*2nkxkzA_A*+e2jk*z 
+k*z2  nkxkzA_B*+ej^+K^z 
-k*.2nkxkzA_A*_ 
+K,*z2^kxkzA_B*_ej{kz~Kz)z 
+k*M-K2x  +  4)B.A*+e^+kz>) 
-«>(-«£  +  K2z)B_B*+e2^z) 
+Kii(~k2x  +  K2z)B_A*_ej{K*-k*)z) 
-k*zii(-4  +  k2z)B_B*_)  . 


(5.34) 
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In  a  solid,  The  individual  intensity  terms  are  not  constant  across  the  layer;  the 
cross-terms  of  the  multiplication  yield  two  position-dependent  terms.  However,  the 
overall  intensity  /  =  Iz  +  Ix  in  either  direction  will  be  constant  in  the  layer. 

Due  to  conservation  of  energy,  the  total  intensity  in  each  material  layer  is  constant 
across  the  system,  and  equal  to  the  outgoing  intensity  in  the  final  layer  In' 

h  —  h  —  In  —  In  ■  (5.35) 

While  the  intensity  is  constant  with  position  across  each  layer,  it  will  vary  as  a 
function  of  frequency  and  incident  angle. 


5.1.2  Model  Verification 


Initial  verification  of  the  model  was  performed  by  matching  the  results  to  closed- 
form  fluid-fluid  and  fluid-solid  interface  solutions  (which  were  presented  in  Chapter 
3),  which  matched  for  all  test  cases,  as  seen  in  Figure  5.3. 

Another  test  of  the  viability  of  this  multi-layer  model  will  be  to  compare  the  model 
results  to  other  models  for  transmission  through  a  simple  plate.  The  transmission 
loss  for  a  plate  calculated  using  Bernoulli-Euler  or  Timoshenko  beam  theory  relies 
on  flexural  wave  motion  in  the  beam;  this  motion  is  simulated  in  the  wave  potential 
model  by  coupling  longitudinal  and  shear  stresses.  The  transmission  loss  across  the 
panel  is: 


T  =  10  log 


1  + 


Z cos (6) 


‘I'PCpl 


(5.36) 


where  the  impedance  Zpi  is  calculated  from  the  plate  flexural  wave  velocity  cpi : 


Zpi  =  jurripi^l  -  ^  .  (5.37) 

For  Bernoulli-Euler  beam  theory,  the  flexural  wave  velocity  can  be  calculated  using 
the  mass  per  unit  area  mpi  and  bending  stiffness  Bpi : 


Imaq.  Comp,  of  Angle  (radians)  Imaq.  Comp,  of  Anqle  (radians) 
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Figure  5.3.  Comparison  of  the  longitudinal  (left)  and  shear  (right) 
transmission  coefficients  for  the  fluid-solid  interface  obtained  using 
the  model  in  Chapter  3  (top  row)  compared  to  that  using  the  multi¬ 
layer  model  (bottom  row).  That  the  plots  along  each  column  match 
is  indicative  of  similar  results  between  both  models. 


Cpl,BE 


(5.38) 


Calculation  of  the  Timoshenko  beam  theory  transmission  loss  uses  a  different  value  of 


cpi ,  which  can  be  found  (as  detailed  can  be  found  in  Graff’s  book  [47])  by  solving  for 


the  positive  root  of  the  frequency  equation  of  the  governing  equations  for  Timoshenko 
beam  motion: 


mpl\cpl,T )  t  V  Ur  J  \Cpl,T )  Gt 


0, 


(5.39) 


where  E,  G,  and  p  are  material  properties  of  the  plate,  /  =  f3/ 12  is  the  bending 
stiffness  of  the  plate,  and  t  is  the  thickness  of  the  plate. 

The  plate  modeled  is  a  3  centimeter  thick  steel  plate,  with  an  acoustic  wave 
incident  at  a  30  degree  angle.  A  plot  showing  the  transmission  loss  variation  with 
frequency  for  the  wave  potential  model,  the  Bernoulli-Euler,  and  Timoshenko  beam 
theories  is  shown  in  Figure  5.4.  Differences  in  the  transmission  loss  are  due  to  the 
assumptions  in  each  beam  theory  that  limit  the  type  of  wave  propagation  in  the  beam. 
The  wave  potential  model  matches  the  prior  theories  well  across  most  frequencies.  The 
model  estimates  the  transmission  loss  minima  near  2,000  Hertz  at  a  frequency  above 
the  Bernoulli-Euler  beam  theory  model  and  below  the  Timoshenko  theory  model. 
Above  the  TL  minima,  the  models  agree  to  within  several  dB,  with  the  wave  potential 
model  more  closely  matching  the  Timoshenko  beam  theory  that  would  allow  for  less- 
constrained  motion  in  the  beam.  Comparison  to  previously-verified  TL  calculations 
shows  that  the  wave  potential  model  is  accurate  for  modeling  acoustical  propagation 
through  materials.  In  addition  to  a  single  layer  of  the  model,  Figure  5.4  shows  two 
additional  lines,  comprising  the  same  panel  modeled  using  two  and  three  sublayers. 
These  lines  lie  atop  the  single-layer  model,  showing  the  robustness  of  the  model. 


5.2  Wavenumber-Frequency  Decomposition  of  Vibration 

Wavenumber-frequency  displays  of  vibration  can  be  useful  for  visualizing  wavetypes 
that  can  propagate  through  a  material.  By  representing  the  vibration  spectrum  in 
the  wavenumber  domain  as  well  as  the  frequency  domain,  properties  of  waves  in  the 
materal  can  be  determined  in  a  way  that  is  not  intuitive  in  the  spatio-frequency 
domain. 
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Figure  5.4.  Transmission  loss  comparison  of  Wave  Potential, 
Bernoulli-Euler  and  Timoshenko  theories,  and  a  split  Wave  Poten¬ 
tial  model  for  2  and  3  sublayers  comprising  the  panel,  for  a  2cm  panel 
with  30  degree  acoustic  incidence. 


Spatio-frequency  vibration  data  V ( x ,  /)  along  a  straight  line  of  length  L,  where 
x  is  the  position  along  the  line  from  0  to  L  and  /  is  the  frequency,  can  be  trans¬ 
formed  into  wavenumber-frequency  data  V ( k ,  /)  by  applying  a  Fourier  transform  to 
the  spatial  dimension.  For  the  case  of  discrete  data  along  equally-spaced  acquisition 
points,  Xi ,  limitations  on  the  Fourier  transform  (as  implemented  using  the  Discrete 
Fourier  Transform)  are  similar  to  those  imposed  in  the  case  of  time-sampled  data 
being  converted  to  frequency.  The  wavenumber  resolution,  A k,  of  a  set  of  discretely 
sampled  points  along  a  length  L  will  be: 

A k  =  27t^— ^  ,  (5.40) 

and  the  maximum  wavenumber  kmax  is  proportional  to  the  spatial  sampling  rate  Ax: 


k 


max 


(5.41) 
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The  resulting  wavenumber-frequency  data  allows  analysis  of  vibration  in  a  do¬ 
main  in  which  wavetypes  are  more  clearly  visible.  High-amplitude  features  that  lie 
along  lines  (either  straight  lines  or  curved)  in  the  wavenumber-frequency  domain  are 
indicative  of  wavetypes  that  are  propagating  in  the  material.  The  phase  speed  of  a 
feature  located  at  (k,u)  is  c  =  uj/k  =  2ir f/k,  and  the  group  speed  is  cg  =  duj/dk.  For 
non-dispersive  waves  which  have  a  speed  invariant  with  frequency,  including  compres- 
sional  and  shear  waves  in  a  material,  the  phase  speed  will  be  constant  and  equal  to 
the  group  speed,  and  the  line  will  be  straight  in  the  wavenumber-frequency  domain. 
For  dispersive  waves,  such  as  flexural  waves  in  a  structure,  the  phase  and  group  speed 
will  vary  with  frequency  and  the  resulting  line  representing  the  wave  will  have  a  vari¬ 
able  slope.  A  wavetype  feature  originating  from  a  f  ^  0  frequency  indicates  that  the 
wave  has  a  cut-on  frequency;  typically,  these  waves  have  geometry  constraint  along 
a  non-propagation  direction  that  must  be  fulfilled  before  the  wave  can  propagate. 
Single  points  of  high  amplitude  are  indicative  of  standing  waves  and  mode  shapes 
in  the  system  (although  mode  shapes  made  from  reflecting  wavetypes  will  show  up 
as  higher-amplitude  sections  of  an  existing  wavetype  feature).  A  schematic  diagram 
showing  these  relationships  on  a  simplified  wavenumber-frequency  plot  is  shown  in 
Figure  5.5. 

The  wavenumber-frequency  decomposition  technique  can  be  applied  to  both  ex¬ 
perimental  measurements  and  theoretical  models.  We  will  use  this  technique  in  dis¬ 
playing  the  results  of  our  model  to  more  easily  distinguish  the  characteristics  in  the 
system. 

5.3  Modeling  of  Wave  Propagation 

Modeling  several  representative  multi-layer  systems  can  simulate  the  effects  of 
coupling  different  materials  together.  Here,  we  build  upon  the  fluid-fluid  and  fluid- 
solid  interaction  by  adding  additional  complexities,  including  a  single  thick  solid  layer, 
the  effects  of  higher-impedance  bounding,  and  the  effects  of  fluid  gaps  between  the 
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Figure  5.5.  Wavenumber-frequency  schematic  showing  phase  speed, 
group  speed,  cut-on  frequency  relationshpis,  and  standing  waves. 


bounding  and  target  layer.  The  model  used  here  allows  for  visualization  of  energy 
within  each  layer,  which  is  useful  for  judging  the  effectiveness  of  difference  incident 
wave  frequencies  and  incident  angles. 

Results  from  the  model  are  a  function  of  the  frequency  /  and  the  incident  angle 
of  wave  propagation,  9.  The  trace  wavenumber  on  the  surface  kx ,  which  is  constant 
throughout  each  layer  of  the  material,  can  be  related  to  the  angle  of  incidence  by: 

kx  =  ksm(9 )  ,  (5-42) 

thus  allowing  for  a  method  to  input  a  surface  wavenumber  into  the  model  using  the 
angle  9.  The  use  of  the  surface  wavenumber  as  a  reference  is  advantageous  because 
vibration  scans  of  the  surface  will  more  closely  correlate  in  that  case.  The  nature  of 
surface  wavenumbers  is  such  that  they  will  vary  from  0  at  normal  incidence  to  the 
incident  region  material  wavenumber  k0  at  90  degrees.  Figure  5.6  shows  equal-angle 
lines  on  a  frequency- angle  plot  and  on  a  wavenumber-frequency  plot.  Each  angle 
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corresponds  to  a  trace  wave  speed  of  c/  8111(6*);  the  decrease  of  trace  wave  speed  with 
increasing  angle  can  clearly  be  seen  in  the  wavenumber- frequency  plot. 


Figure  5.6.  Representation  of  incident  angles  of  propagation  in 
wavenumber-frequency  plots.  Lines  of  equal  angle  are  shown  in  both 
the  frequency-angle  plot  (left)  and  frequency- wavenumber  plot  (right). 


Because  the  surface  wavenumber  is  a  function  of  frequency,  the  variables  in  the 
wavenumber  frequency  plot  are  not  independent.  This  leads  to  compression  of  some 
regions  in  the  plot,  particularly  at  lower  wavenumbers,  where  the  entire  angle  is 
compressed  into  a  much  smaller  section  than  the  equivalent  angle  range  at  higher 
frequencies. 

Several  materials  will  be  simulated  in  the  model  which  are  close  analogues  to  the 
materials  of  interest  in  our  testing  and  research.  These  properties  were  arrived  at 
through  a  combination  of  existing  data  and  wave  propagation  testing  on  the  material 
(as  detailed  in  Section  6),  and  are  adjusted  to  yield  properties  that  can  be  clearly 
distinguished  in  analysis  when  coupled  to  other  materials. 

Our  primary  concern  is  a  low-wave  speed  material  that  approximates  an  energetic 
material.  The  material  is  question  is  a  composite  polymer  of  ammonium  chloride  in  a 
hydroxyl-terminated  polybutadiene  (HTPB)  resin  binder,  with  the  ammonium  chlo- 
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ride  consisting  of  75%  of  the  polymer  by  weight.  Further  details  of  the  material  will 
be  given  in  Chapter  6;  for  our  purposes  now,  we  approximate  the  material  properties 
as  longitudinal  and  shear  wave  speeds  of  700  and  500  meters  per  second,  respectively, 
and  a  density  of  1,200  kilograms  per  cubic  meter.  Because  our  primary  goal  is  to 
excite  this  energetic  material,  results  of  the  system  will  be  shown  in  terms  of  energy 
deliverance  (the  sum  of  both  longitudinal  and  shear  intensity)  into  the  target  layer 
of  surrogate  material. 

Polycarbonate  is  used  as  a  bounding  material,  with  a  longitudinal  wave  speed  of 
1,300  meters  per  second,  a  shear  wave  speed  of  900  meters  per  second,  and  a  density 
of  1,000  kilograms  per  cubic  meter.  The  fluid  in  the  model  is  air,  with  a  wave  speed  of 
343  meters  per  second  and  a  density  of  1.2  kilograms  per  cubic  meter.  The  properties 
used  in  all  models  are  shown  in  Table  5.1. 

Table  5.1.  Table  of  material  properties  used  in  multi-layer  wave 
propagation  model. 


Material 

Density 

P 

kg/m3 

Compressional 

wave  speed 

c 

m/s 

Shear 

wave  speed 

b 

m/s 

Air 

1.2 

343 

N/A 

Surrogate 

1,200 

700 

500 

Polycarbonate 

1,000 

1,300 

900 

To  limit  the  effect  of  resonances  in  the  layers,  a  1%  damping  coefficient  was  added 
to  the  speed  of  sound  of  all  enclosed  layers. 

The  initial  model  we  will  consider  is  a  fluid  layer  of  10  centimeters  bounded  by  air 
on  both  sides,  with  the  density  and  (longitudinal)  wave  speed  of  the  fluid  equal  that 
of  the  energetic  material.  As  in  Chapter  3,  this  will  serve  as  a  simplified  model  of 
wave  propagation  to  allow  for  ease  of  analysis.  This  same  geometry  with  a  solid  layer 
of  energetic  material  will  also  be  considered.  In  the  next  model,  a  thin,  2  millimeter 
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layer  of  polycarbonate  will  border  the  energetic  material  to  provide  a  high-impedance 
material  barrier  to  wave  transmission.  Lastly,  the  effects  of  having  a  fluid  gap  in  the 
system  between  the  barrier  and  the  energetic  material  will  be  explored  by  inserting  a 
1  centimeter  gap  of  air  into  the  system.  A  schematic  of  all  the  multi-layer  geometries 
are  shown  in  Figure  5.7. 

Unless  noted,  the  intensity  transmission  loss  is  plotted  in  the  figures  below: 

TL  =  mog(^^j  ,  (5.43) 

where  z  is  typically  the  center  of  the  i-th  layer  of  interest  (which  is  noted  in  each 
individual  analysis). 

The  wavenumber-frequency  characteristics  of  the  intensity  transmission  loss  at 
the  center  of  the  single  fluid  layer  is  shown  in  Figure  5.8.  The  region  that  comprises 
subsonic  transmission  in  air  (below  343  meters  per  second  of  wave  speed)  has  been 
excluded  because  the  resulting  incident  wave  generated  in  this  region  would  not  prop¬ 
agate,  and  the  reference  intensity  is  too  low  to  provide  a  physically  relevant  measure 
of  transmision  loss  in  the  system. 

Between  the  wave  speed  of  the  incident  air  and  the  lowest  wave  speed  in  the  mate¬ 
rial  of  interest  (in  this  case,  the  fluid  wave  speed  of  700  meters  per  second),  all  waves 
into  the  material  are  subsonic;  the  region  of  high  transmission  loss  here  is  indicative 
of  very  little  energy  propagating  into  the  material.  Although  the  temperature  plot 
of  the  transmission  loss  is  cut  off  at  100  dB  to  show  detail  in  the  supersonic  region, 
the  actual  transmission  loss  in  this  region  is  around  200  dB,  and  will  be  of  no  use  in 
energy  transmission  into  the  material. 

At  k  —  0,  which  corresponds  to  normal  incidence  of  sound,  there  are  modes  at  in¬ 
teger  multiples  of  3,500  Hertz,  corresponding  to  the  half-wavelength  distance  of  sound 
in  the  energetic  material  /  =  racj/Lj.  With  increasing  wavenumber,  the  frequency 
of  these  features  increase  in  frequency  to  account  for  the  longer  distance  imposed  on 
the  half- wavelength.  They  eventually  become  waves  moving  almost  entirely  down  the 
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Polycarbonate  Walls 
2  mm  2  mm 


Figure  5.7.  Geometries  of  modeled  multi-layer  systems.  The  top 
system  consists  of  10  centimeters  of  a  surrogate  layer  (and  will  also 
be  used  for  the  single  fluid  case).  The  second  system  is  the  same 
surrogate  layer  enclosed  by  thin  polycarbonate  walls,  to  simulate  a 
container.  The  lowest  geometry  is  the  same  surrogate  and  walls,  only 
with  a  1  centimeter  air  gap  between  the  surrogate  and  the  polycar¬ 
bonate  wall. 


waveguide  of  the  center  fluid  layer,  and  the  group  speed  of  waves  from  this  motion 
will  approach  the  speed  of  longitudinal  wave  propagation  in  the  material. 
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Figure  5.8.  Intensity  in  middle  layer  of  the  fluid  layer  system,  with 
the  wave  speed  of  the  fluid  (700  meters  per  second)  and  the  wave 
speed  in  air  (343  meters  per  second)  indicated. 


In  the  results  for  a  solid  layer  of  material,  shown  in  Figure  5.9,  the  addition  of 
shear  waves  to  the  system  adds  several  features.  The  high-transmission  loss  region 
seen  in  the  fluid  layer  modeled  is  not  seen  in  this  model;  while  there  ought  to  be  a 
high-transmission  loss  region  below  the  shear  speed  in  the  material,  that  region  is 
dominated  by  a  flexural  mode  which  decreases  transmission  loss.  The  effect  of  the 
shear  wave  on  energy  transmitted  into  the  material  is  visible  because  of  the  features  at 
integer  multiples  of  2,500  Hertz,  which  is  the  half  wavelength  frequency  of  materials 
traveling  at  500  meters  per  second  in  a  10  centimeter  gap. 

Above  normal  incidence,  the  complex  behavior  introduced  into  the  system  by  the 
addition  is  due  to  wave  coupling  phenomenon,  wherein  the  shear  waves  that  were  not 
present  in  the  fluid  layer  couple  with  the  existing  longitudinal  waves  to  create  Lamb 
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Figure  5.9.  Intensity  in  the  middle  layer  of  the  solid  layer  system, 
with  relevent  features  and  wave  speeds  labeled. 


waves  propagating  within  the  solid  layer.  This  causes  behavior  such  as  negative  group 
speed  seen  in  the  curve  beginning  at  3,500  Hz  (Feature  1).  Auld  [22]  goes  into  this 
phenomenom  more  fully. 

The  addition  of  the  polycarbonate  bounding  layer  serves  to  lower  the  cut-on  fre¬ 
quency  and  the  frequencies  for  the  equivalent  wavetypes  in  the  material,  as  shown  in 
Figure  5.10.  This  is  more  evident  at  low  frequencies,  where  the  mode  previously  found 
at  3.5  kilohertz  is  now  at  2.8  KHz.  While  the  polycarbonate  is  a  higher- impedance 
material,  the  layer  is  too  thin  to  contribute  additional  stiffness,  while  increasing  the 
length  of  the  overall  system  in  which  waves  can  propagate. 

The  addition  of  an  air  gap  between  the  two  materials  greatly  increases  transmission 
loss  across  all  angles  and  wavenumbers,  as  seen  in  Figure  5.11.  This  shows  the  effects 
of  decoupling  the  polycarbonate  layers  from  the  energetic;  when  the  polycarbonate 
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Figure  5.10.  Intensity  in  the  center  of  the  middle  (surrogate)  layer 
of  the  plastic-bounded  surrogate  system. 


was  coupled  to  the  solid  layer  as  shown  in  Figure  5.10,  the  resulting  system  had 
little  effect  on  the  transmission  loss.  However,  decoupling  the  layers  greatly  increases 
the  transmission  loss,  as  the  polycarbonate  layer  induces  a  significant  loss  before 
transmission  into  the  energetic  layer. 

5.3.1  The  Effect  of  Evanescent  Incident  Waves 

Evanescent  waves  were  shown  to  be  effective  for  overcoming  subsonic  transmission 
constraints  in  Chapter  3;  their  effect  in  multi-layer  systems  will  allow  for  additional 
avenues  of  energy  transmission. 

In  the  model  developed  here,  evanescent  waves  are  represented  by  an  adding  an 
imaginary  component  to  the  trace  surface  wavenumber.  Such  a  wavenumber  will 
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Figure  5.11.  Intensity  in  the  middle  layer  of  middle  (surrogate)  layer 
of  the  plastic-bounded  system  with  an  air  gap  between  the  plastic  and 
surrogate. 


relate  to  a  complex  incident  angle  of  propagation  using  our  previous  formulation  in 
Equation  5.42. 

The  evanescent  waves  will  not  necessarily  increase  transmission  into  the  system; 
rather,  they  will  allow  for  additional  transmission  in  subsonic  (below  the  material 
wave  speed)  regions  of  the  materials.  For  example,  the  intensity  transmission  loss 
due  to  a  0.1  complex  component  added  to  the  single  solid  layer  shown  in  Figure  5.9 
is  shown  on  the  left  side  of  Figure  5.12.  The  scale  of  this  plot  has  been  left  at  the 
same  scale  of  the  previous  plot;  the  right  plot  shows  the  reduction  in  transmission 
loss  between  the  transmission  loss  of  the  evanescent  wave,  TLe  compared  to  the 
transmission  loss  of  the  classical  plane  wave,  TLp,  TLp  —  TLe.  In  such  a  case,  a 
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positive  value  of  TL  reduction  indicates  an  increase  in  energy  transmission  due  to 
evanescent  wave  effects. 

The  use  of  evanescent  waves  causes  increased  energy  transmission  into  the  ma¬ 
terial,  particularly  at  lower  wave  speeds,  up  to  20  dB.  The  evanescent  wave  also 
broadens  the  region  of  low  transmission  loss  (Feature  2),  increasing  the  width  of  the 
flexural  mode  components  to  lower  the  transmission  loss  in  that  region.  There  is  an 
additional  narrow  high-transmission  loss  region  that  is  below  a  surface  wavenumber 
of  100  radians  per  meter  (Feature  1);  the  cause  of  this  is  unknown. 
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Figure  5.12.  Intensity  transmission  loss  in  the  center  layer  (left)  and 
transmission  loss  improvement  due  to  the  evanescent  wave  (right) 
for  the  single-layer  surrogate  layer  system  with  0.1  radian  evanescent 
component. 


Evanescent  waves  ought  to  provide  a  method  for  propagating  energy  across  a 
higher  wave  speed  barrier  region,  particularly  when  the  energetic  behind  the  barrier 
has  a  lower  wave  speed.  For  example,  a  thick  polycarbonate  layer  bonded  to  the 
surrogate  material  will  provide  significant  energy  loss  in  large  regions  of  transmission. 
Such  a  system,  with  a  5  centimeter  layer  of  polycarbonate  in  front  of  10  centimeter 
layer  of  surrogate  material,  will  show  the  contribution  of  evanescent  waves.  The  in- 
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tensity  transmission  loss  into  the  center  of  the  polycarbonate  and  surrogate  layers 
for  non-evanescent  wave  transmission  is  shown  in  Figure  5.13.  Intensity  transmission 
in  the  system  is  dominated  by  the  polycarbonate  layer,  causing  subsonic  transmis¬ 
sion  into  the  material  below  700  meters  per  second.  The  addition  of  an  evanescent 
component,  shown  in  Figure  5.14,  shows  regions  of  improved  and  decreased  intensity 
transmission  loss.  Figure  5.15  shows  the  TL  reduction  in  the  center  layer  due  to  the 
evanescent  component.  At  higher  wave  speeds,  there  is  a  decrease  in  intensity  trans¬ 
mission  loss  (Feature  1);  this  may  be  due  to  the  increased  decay  of  wave  propagation 
limiting  total  intensity  to  the  layer.  At  lower  wave  speeds  (Feature  2),  there  is  a 
slight  increase  in  energy  transmission  (about  10  dB);  this  could  be  a  region  where 
the  additional  decay  in  the  system  is  still  greater  than  the  effects  of  classical  wave 
propagation,  allowing  a  contribution  from  the  evanescent  wave. 
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Figure  5.13.  Intensity  transmission  loss  for  the  polycarbonate- 
shielded  surrogate  system,  consisting  of  polycarbonate  in  front  of  10 
centimeters  of  surrogate.  Intensity  transmission  loss  is  shown  in  the 
polycarbonate  layer  (left)  and  the  surrogate  layer  (right). 
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Figure  5.14.  Intensity  transmission  loss  with  evanescent  component 
for  the  polycarbonate-shielded  surrogate  system,  consisting  of  poly¬ 
carbonate  in  front  of  10  centimeters  of  surrogate.  Intensity  transmis¬ 
sion  loss  is  shown  in  the  polycarbonate  layer  (left)  and  the  surrogate 
layer  (right). 
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Figure  5.15.  Intensity  transmission  reduction  in  the  surrogate  layer 
for  the  polycarbonate-shielded  surrogate  system. 
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5.4  Conclusions 

Wave  propagation  in  multi-layer  materials  is  governed  by  the  coupled  nature  of 
the  component  layers.  While  previous  models  have  developed  transfer  matrix  and 
system  of  equation  formulations  for  wave  propagation  through  a  multi-layer  system, 
using  wave  potentials  to  calculate  wave  propagation  has  several  advantages,  such  as 
allowing  for  evanescent  incident  waves  and  calculation  of  intensity  in  the  system. 

The  wave  potential  model  is  generated  for  modeling  wave  propagation  through 
each  component  layer  of  material,  and  equating  the  material  states  at  each  interface 
to  solve  for  the  appopriate  wave  properties.  By  replacing  the  transfer  matrix-adapted 
method  of  wave  propagation  with  a  wave  pontential  method,  the  material  states  can 
be  seen  and  intensity  in  each  layer  can  be  calculated.  The  model  was  verified  by 
matching  the  results  to  existing  closed-form  solutions  of  panel  vibration. 

Several  multi-layer  materials  were  modeled,  beginning  with  a  single  fluid  or  solid 
layer.  A  single  layer  exhibits  both  modes  that  are  integer  wavelengths  of  the  layer’s 
thickness,  as  well  as  coupled  modes  that  are  flexural  or  dilatational  motion  in  the 
solid.  When  bounded  by  a  higher-impedance  material  with  low  flexural  stiffness,  the 
materials  couple  together  and  produce  a  lower  flexural  stiffness  multi-layer  material. 
The  removal  of  the  coupling  between  the  layers  by  adding  a  small  fluid  layer  between 
each  solid  produces  a  much  higher  energy  loss  across  the  sample. 

The  model  allows  for  evanescent  incident  waves;  these  effects  are  seen  in  the 
subsonic  region  of  energy  transmission,  but  have  little  or  no  effect  in  the  classic 
regions.  An  evanescent  waves  can  be  used  to  propagate  energy  into  the  subsonic 
domain  of  a  material,  which  can  be  useful  to  overcome  a  high-impedance  barrier  before 
a  low-impedance  target  material.  In  cases  where  a  particular  frequency  is  desired 
that  may  impose  limitations  on  classical  wave  propagation,  the  use  of  evanescent 
angles  can  add  another  parameter  to  overcome  those  limitations.  However,  the  model 
constructed  shows  many  regions  of  classical  wave  propagation  that  are  as  effective  or 
better  to  excite  material  layers. 
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Because  the  evanescent  waves  offer  improvement  in  sound  transmission  in  the  sub¬ 
sonic  region,  they  pale  in  comparison  to  intensity  transmission  in  supersonic  regions. 
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CHAPTER  6.  EXPERIMENTAL  MEASUREMENTS  OF  WAVE 

PROPAGATION 

The  ability  to  quantify  wave  propagation  characteristics  in  materials  and  structures 
experimentally  depends  significantly  on  conditions  inherent  to  the  testing  setup.  The 
resulting  wave  speeds  found  in  testing  can  be  used  as  inputs  into  the  multi-layer  model 
described  in  Chapter  5,  as  well  as  verification  for  that  same  model  (particularly  for 
geometry-based  considerations  such  as  flexural  waves).  Direct  measurement  of  wave 
propagation  is  particularly  important  for  inhomogeneous  materials  such  as  foams  or 
aggregate  composites,  which  may  have  wave  properties  that  differ  significantly  from 
what  the  material  properties  would  suggest.  Furthermore,  measurement  of  propaga¬ 
tion  characteristics  will  provide  a  better  estimate  of  properties  such  as  the  speed  of 
wave  propagation,  as  well  a  geometry-specific  properties  such  as  cut-on  frequencies 
or  flexural  wave  speeds. 

Vibration  visualized  in  the  wavenumber-frequency  domain  has  distinct  features 
that  correspond  to  properties  such  as  wave  speeds,  modes,  and  cut-on  frequencies  in  a 
component.  However,  vibration  measured  during  testing  will  be  subject  to  conditions 
of  the  test,  in  particular  excitation  position  and  mounting  conditions.  Poor  mounting 
conditions  will  mask  the  resulting  waves  in  the  system;  while  free-free  simulated 
mounting  commonly  used  in  vibration  classification  of  materials  provide  an  ideal 
condition,  they  are  insufficient  in  cases  where  the  material  needs  support,  and  the 
mounting  must  be  considered. 

Another  factor  to  be  considered  in  experimental  measurements  is  the  effect  of  ex¬ 
citation  locations,  which  may  reveal  different  wavetypes  in  the  system.  For  example, 
excitation  at  the  end  of  a  beam  may  excite  longitudinal  waves  down  the  bar  at  the 
expense  of  transverse  excitation.  Care  should  be  taken  in  designing  an  experiment  to 
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assess  the  type  of  input  and  the  effects  it  may  have  on  the  material.  While  the  excita¬ 
tion  type  may  be  useful  for  illustrating  different  wavetypes,  determining  their  relative 
susceptibility  to  excitation  should  only  be  performed  with  careful  consideration  of  the 
system. 

In  this  chapter,  the  effects  of  mounting  conditions  and  excitation  location  are 
considered  in  the  wave  propagation  classification  of  an  unknown  material. 

6.1  Test  Methodology 

The  motivation  of  this  research  is  determining  the  wave  propagation  characteris¬ 
tics  of  a  surrogate  explosive  material.  Other  methods  of  classification  are  being  pur¬ 
sued  simultaneously  by  Paripovic  and  Davies  [48];  however,  direct  measurement  of 
wave  properties  will  be  useful  for  modeling  purposes.  The  material  of  interest  is  com¬ 
posite  polymer  composed  of  an  ammonium  chloride  crystal  in  a  hydroxyl-terminated 
polybutadiene  (HTPB)  resin  binder,  mixed  together  into  a  homogeneous  mixture  and 
poured  into  molds  where  they  hardened  to  create  a  solid  material.  The  surrogate  is 
classified  by  the  percentage  of  the  crystal  volume  fraction,  with  50%  and  75%  samples 
being  used  in  our  testing.  The  material  was  designed  this  way  to  simulate  the  dy¬ 
namic  properties  of  a  commonly-used  propellant  while  avoiding  rapid  energy  release 
should  the  material  be  over-excited.  In  the  propellant,  the  ammmonium  chloride 
crystals  would  be  replaced  by  ammonium  perchloride. 

The  HTPB  samples  are  cast  into  a  mold  1.27  centimeters  (0.5  inch)  thick  by  2.54 
centimeters  (1  inch)  wide  and  121.9  centimeters  (48  inches)  long.  The  properties  of 
the  HTPB  surrogate  vary  significantly  between  the  50%  and  75%  samples.  While  the 
50%  samples  are  flexible  and  deform  under  their  own  weight,  the  75%  samples  are 
more  rigid  but  display  significant  granularity  that  is  susceptible  to  breaking  under 
deformation.  Neither  sample  is  rigid  enough  to  stay  intact  if  subject  to  its  own  weight, 
and  cannot  be  machined  or  tapped  for  attachment  of  measuring  devices.  Figure  6.1 
is  a  photo  showing  the  textural  differences  between  the  two  samples. 
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Figure  6.1.  Surrogate  samples  tested.  The  50%  sample  is  on  top, 
consisting  of  lower  stiffness  and  a  more  elastic  structure.  The  75% 
sample  is  significantly  more  rigid,  and  would  fracture  under  substan¬ 
tial  deformation. 


Because  the  surrogate  samples  are  fragile  and  cannot  support  their  own  weight, 
mounting  conditions  are  limited,  especially  those  which  provide  a  minimal  support. 
In  order  to  minimize  the  effects  of  the  support,  the  samples  were  placed  on  Thinsnlate 
layer  that  rested  on  a  Quash  foam  board.  A  polycarbonate  backing  was  also  provided 
for  comparison. 

In  addition  to  the  HTPB  surrogate  material,  polycarbonate  was  used  as  well  to 
provide  both  a  low-damping,  thoroughly-quantified  material.  The  polycarbonate  is 
extruded  as  a  thickness  of  1.27  centimeters  (0.5  inch)  and  cut  to  a  width  of  2.54 
centimeters  (1  inch)  and  a  length  of  121.9  centimeters  (48  inches).  The  clear  sides  of 
polycarbonate  were  sanded  to  provide  an  opaque  surface  that  could  be  measured  by 
the  laser  vibrometer. 
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The  density  of  polycarbonate  is  1,200  kilograms  per  cubic  meter  (kg/m3).  The 
50%  surrogate  material  has  a  density  of  1,156  kg/m3,  and  the  75%  has  a  density 
of  1,085  kg/m3.  While  the  ammonium  chloride  material  has  a  higher  density  than 
the  resin  binder,  manufacturing  of  the  75%  samples  produces  air  gaps  that  serve  to 
reduce  the  density  of  the  material.  All  of  the  foams  have  significantly  lower  densities 
than  the  test  materials;  the  polyamide  foam  used  in  the  minimal  foam  backing  tests 
of  polycarbonate  had  a  density  of  6.4  kg/m3,  while  the  Quash  had  a  density  of  32 
kg/m3.  The  volumetric  density  of  the  Thinsulate  could  not  be  found;  however,  the 
sheet  material  weighs  150  kg/m2  for  a  sheet  approximately  6  millimeters  thick  when 
uncompressed,  which  would  yield  a  density  of  9  kg/m3.  The  density  of  the  tested 
materials  is  at  least  30  times  greater  than  the  backing  material,  and  typically  more. 

The  polycarbonate  samples  were  tested  in  three  mounting  conditions.  The  first 
was  a  hung  condition,  wherein  the  polycarbonate  rod  was  suspended  from  a  hook  by 
a  string  drilled  into  the  far  edge  of  the  sample.  This  represents  a  minimal  boundary 
condition,  and  is  commonly  used  in  modal  testing  of  structures.  The  polycarbonte 
was  then  supported  using  the  minimal  foam  support,  as  well  as  a  more  thorough  foam 
backing  consisting  of  the  Thinsulate,  the  same  method  used  for  the  surrogate  sample. 

The  excitation  locations  used  in  testing  were  a  side  excitation,  wherein  the  actu¬ 
ator  was  mounted  within  an  inch  from  the  end  of  the  largest  side  (on  the  width-by- 
length  plane).  The  other  excitation  location  was  on  the  end  (on  the  width-by-height 
plane).  The  actuator-transducer  stack  was  attached  to  the  polycarbonate  through  a 
screw  stud  that  was  tapped  into  the  polycarbonate.  To  attach  the  stack  to  the  HTPB 
surrogate  samples,  an  accelerometer  mount  was  glued  to  the  sample  surface  by  using 
epoxy;  these  mounts  have  a  screw  stud  that  allows  the  stack  to  be  attached. 

The  samples  were  excited  using  a  PCB  712A02  stack  actuator  equipped  with  a 
100  gram  add-on  mass.  The  actuator  has  an  excitation  range  of  150  to  5000  Hertz, 
and  was  driven  with  a  white- noise  signal  generated  by  a  General  Radio  1381  Random- 
Noise  Generator,  running  at  a  brandwidth  of  5  KHz.  The  signal  was  amplifed  with  a 
QSC  1100  audio  amplifier. 
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Figure  6.2.  Schematic  of  experimental  measurements  of  polycar¬ 
bonate  beam.  Two  excitation  locations  were  used,  with  the  ’’side” 
excitation  being  along  the  same  plane  as  the  measurement,  and  the 
’’end”  excitation  being  at  the  far  end  of  the  polycarbonate  bar.  For 
mounting,  the  beam  was  hung  (suspended  from  the  end  above  a  fixed 
point),  given  a  foam  backing  that  covered  the  entirety  of  the  surface 
opposite  the  measurement,  and  positioned  with  minimal  support  that 
allowed  for  reduced  contact  with  the  sample. 


Force  into  the  sample  was  measured  using  a  PCB  208A02  force  transducer,  which 
is  an  ICP-powered  transducer  with  a  nominal  output  of  50  millivolts  per  Newton. 
This  transducer  was  positioned  between  the  actuator  and  the  sample. 

Vibration  measurement  along  the  sample  was  performed  using  a  Polytec  PSV-400 
scanning  laser  vibrometer.  The  vibrometer  measured  a  4  KHz  bandwidth  with  1,600 
points  across  that  frequency  range  for  a  frequency  resolution  of  2.5  Hz.  The  number 
of  points  tested  varied  during  particular  tests;  as  least  200  points  where  measured 
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across  every  sample.  Testing  spanned  over  most  of  the  bar,  with  a  scan  length  of 
approximately  1.15  meters  in  length. 

All  data  analysis  was  performed  using  the  data  acquisition  system  of  the  Polytec 
PSV-400;  both  the  vibration  signal  and  the  force  transducer  signal  were  windowed 
using  a  Hann  window  and  recorded  for  50  averages  with  50  percent  overlap.  The  HI 
estimator  was  calculated  between  the  vibration  and  force  transducer  signal,  and  is 
the  displayed  signal  on  all  results.  A  list  of  all  the  equipment  used  in  testing  is  shown 
in  Table  6.1. 

Wave  propagation  measurement  consists  of  exciting  vibration  in  the  sample  using 
the  actuator,  then  measuring  the  resulting  vibration  at  number  of  points  along  the 
length  of  the  sample.  Conversion  of  the  frequency-position  data  into  wavenumber- 
frequency  data  enables  visualization  of  wavetypes  in  the  system.  Wavenumber  anal¬ 
ysis  was  performed  using  the  method  specified  in  Section  5.2.  A  table  showing  exper¬ 
imental  parameters  and  the  resulting  frequency  and  wavenumber  resolution  is  shown 
in  Table  6.2. 

Use  of  the  Discrete  Fourier  Transform  results  in  a  two-sided  wavenumber-frequency 
plot,  with  information  in  both  the  negative  and  positive  wavenumber  values.  This 
information  is  indicative  of  waves  traveling  in  different  directions  in  the  system;  the 
direction  is  relative  to  the  numbering  of  points  in  the  system,  and  in  all  cases  points 
are  numbered  with  low  numbers  being  closer  to  the  excitation  location,  and  waves 
traveling  away  from  the  sensor  have  negative  wavenumbers.  Because  waves  traveling 
towards  the  excitation  source  can  only  be  formed  by  reflection  in  the  material,  the 
difference  between  the  features  on  the  positive  and  negative  sides  of  the  wavenumber 
spectrum  is  due  to  reflection  characteristics  in  the  system. 
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Table  6.1.  Equipment  used  in  wave  propagation  testing. 


Equipment  Type 

Manufacturer 

Model  Number 

Serial  Number 

Stack  Actuator 

PCB 

712A02 

666 

White  Noise  Source 

General  Radio 

1381 

41259 

Amplifier 

QSC 

1100 

109298995 

Force  Transducer 

PCB 

208A02 

9273 

Scanning  Vibrometer 

Polytec 

PSV-400 

0111421 

Data  Acquisition 

Polytec 

PSV-400 

0111891-01 

Table  6.2.  Parameters  for  sample  testing. 


Parameter 

Quantity 

Units 

Measurement  Distance 

1.15 

meters 

Spatial  Resolution 

5.6 

millimeters 

Wavenumber  Resolution 

5.3 

rad/m 

Sampling  Frequency 

4,000 

Hz 

Frequency  Resolution 

2.5 

Hz 

Number  of  Averages 

50 

N/A 

Overlap 

50 

percent 

6.2  Results 

6.2.1  Polycarbonate  Tests 

The  results  of  the  tests  on  the  polycarbonate  samples  are  shown  in  Figure  6.3. 
In  the  foam  mounting  condition,  the  vibration  spectrum  of  the  side  excitation  is 
dominated  by  a  dispersive  wave  with  a  phase  speed  between  114  meters  per  second 
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at  500  Hz  and  332  meters  per  second  at  3,500  Hz.  This  is  most  likely  a  flexural  wave 
with  deformation  normal  to  the  direction  of  measurement  and  excitation. 

The  end  excitation  shows  the  flexural  wave,  but  has  an  equally  prominent  non- 
dispersive  wave  of  about  1,700  meters  per  second.  This  would  most  likely  be  the 
longitudinal  wave  speed  in  the  beam;  from  the  material  properties,  polycarbonate  is 
estimated  to  have  a  bulk  longitudinal  wave  speed  of  1,696  meters  per  second.  While 
the  flexural  wave  is  seen  in  the  end  excitation,  another  slightly  faster  dispersive  wave 
can  also  be  observed.  This  is  the  transverse  flexural  wave  causing  motion  in  the 
normal  direction.  The  wave  speed  of  a  flexural  wave  is  proportional  to  the  fourth- 
root  of  the  bending  moment;  the  approximate  ratio  between  the  two  speeds,  across  all 
frequencies,  is  1.7,  which  is  nearly  equal  to  the  ratio  of  normal  to  transverse  bending 
moment,  81/4  =  1.68. 

Both  the  foam  backed  and  minimal  foam  support  mounting  conditions  show  the 
same  trend  with  regards  to  excitation  as  the  hung  mounting,  wherein  the  side  excita¬ 
tion  exhibits  a  strong  flexural  wave  and  the  end  excitation  shows  multiple  wavetypes. 
In  contrast  to  the  hung  mounting  condition,  though,  both  foam  mounting  conditions 
exhibit  the  transverse  flexural  wave  more  prominently.  This  would  indicate  that  the 
normal  flexural  wave  is  limited  by  the  mounting  conditions,  while  the  transverse  flex¬ 
ural  motion  is  not  subject  to  the  same  constraint.  The  additional  prominence  of 
the  transverse  wave  in  the  minimal  foam  support  condition  is  due  to  the  decreased 
restriction  in  motion  of  the  minimal  foam  support  compared  to  the  foam  backing. 
The  structure  of  the  non-dispersive  wave  shown  in  the  minimal  foam  supports  with 
end  excitation  suggests  that  there  are  two  waves  with  this  high  wave  speed;  however, 
it  is  unknown  what  this  wavetype  could  be,  and  could  be  additional  energy  of  the 
longitudinal  wavetype. 

The  tests  above  show  the  results  of  changing  the  excitation  location  and  mounting 
type  on  a  polycarbonate  beam.  The  hung  mounting  provides  the  best  excitation  of  the 
wavetypes  typically  normal  to  the  direction  of  measurements,  but  the  foam  mountings 
reveal  additional  wavetypes  such  as  the  transverse  flexural  wave.  Similarly,  the  end 
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Figure  6.3.  Polycarbonate  vibration  for  the  hung  (top  row),  foam 
backing  (middle  row)  and  minimal  foam  support  (bottom  row)  con¬ 
ditions,  with  side  (left  plots)  and  end  (right  plots)  shown.  Significant 
features  have  been  labeled  on  the  plot. 
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excitation  shows  several  low  level  wavetypes,  while  a  side  excitation  is  dominated  by 
a  single  wave. 

6.2.2  Surrogate  Tests 

The  surrogate  sample  differs  from  the  polycarbonate  in  that  the  material  proper¬ 
ties  were  not  known,  making  estimation  and  suggestion  of  what  the  visualize  wavetypes 
might  be  attributed  to  in  the  material  difficult.  Therefore,  the  polycarbonate  results 
will  provide  a  template  for  which  to  interpret  the  features  found  in  the  surrogate 
tests. 

Figure  6.4  shows  the  vibration  spectra  for  the  50%  surrogate  sample  tests.  The 
sample  exhibits  a  one-sided  wave  spectrum;  this  is  indicative  of  a  lack  of  reflection 
in  the  material  due  to  damping.  The  effects  of  the  damping  can  also  be  seen  in  the 
frequency  response;  there  is  substantially  less  energy  above  2,500  Hertz.  The  lack 
of  energy  at  higher  frequencies  can  be  seen  in  Figure  6.5,  which  shows  the  position- 
dependent  frequency  response.  The  higher-frequency  components  decay  quickly,  and 
therefore  contribute  less  energy  in  the  wavenumber  domain. 

The  50%  surrogate  test  results  show  three  waves  of  interest.  Two  dispersive  waves 
in  the  same  normal-transverse  pair  seen  in  the  polycarbonate  have  maximum  speeds 
(before  damping  becomes  significant)  of  95  and  113  meters  per  second,  are  visible  in 
the  end  excitation,  while  only  the  slower  wave  can  be  seen  in  the  side  excitation.  The 
non-dispersive  wave  has  a  wave  speed  of  244  meters  per  second,  which  is  significantly 
slower  than  the  sound  speed  in  air.  Therefore,  the  50%  surrogate  material  would  have 
no  subsonic  range  of  acoustical  transmission. 

When  the  surrogate  sample  is  backed  with  the  higher-impedance  material,  such 
as  in  the  results  for  the  test  on  the  polycarbonate-backed  surrogate,  the  resulting 
vibration  information  is  indiscernable;  coupling  with  the  polycarbonate  obscures  all 
pertinent  wave  information.  This  would  suggest  that  the  foam  mounting  is  sufficiently 
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Figure  6.4.  50  %  surrogate  vibration  for  foam  backing  (top)  and 

polycarbonate  backing  (bottom)  conditions,  with  side  (left  plots)  and 
end  (right  plots)  shown.  Significant  features  have  been  labeled  on  the 
plot. 


low  impedence  to  allow  for  wave  propagation  in  the  materials,  but  the  polycarbonate 
material  is  too  stiff  and  restrains  motion. 

For  the  results  of  the  75%  surrogate  sample  tests,  shown  in  Figure  6.6,  the  foam 
mounting  conditions  have  much  less  prevalent  damping  than  that  which  produced  a 
one-sided  wavenumber-frequency  profile  for  the  50%  surrogate  sample.  While  the 
forward-going  (left)  side  of  the  plot  is  slightly  higher  in  amplitude,  there  is  still 
significant  energy  traveling  in  the  opposite  direction  through  the  bar.  Two  waves 
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Figure  6.5.  Vibration  in  the  position-frequency  domain  for  the  50% 
surrogate  sample,  comparing  the  foam  mounting  condition  for  side 
(left)  and  end  (right)  excitation. 


are  evident  in  testing:  a  dispersive  wave  with  a  speed  of  about  400  meters  per  second 
above  1000  Hertz,  and  a  non-dispersive  wave  with  a  speed  of  about  900  meters  per 
second.  While  the  transverse  flexural  wave  could  possibly  be  seen  as  the  higher 
component  of  the  main  dispersive  wave,  it  is  not  discernable  enough  from  the  primary 
flexural  wave  to  ascertain  its  properties. 

The  polycarbonate-backed  tests  of  the  75%  sample  have  more  wave  components 
than  the  foam  mounting  condition  and  the  primary  wavetypes  are  harder  to  discern; 
however,  wave  components  can  more  clearly  be  seen  than  in  the  equivalent  test  of 
the  50%  polycarbonate-mounted  sample.  This  suggests  that  the  additional  stiffness 
in  the  75%  sample  makes  it  more  suitable  for  higher- impedance  mounting. 

6.3  Conclusions 

Testing  of  material  samples  yielded  different  results  depending  on  the  excitation 
location  and  mounting  of  the  specimen.  Certain  excitation  effects,  such  as  the  promi- 
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Figure  6.6.  75  %  surrogate  vibration  for  foam  backing  (top)  and 

polycarbonate  backing  (bottom)  conditions,  with  side  (left  plots)  and 
end  (right  plots)  shown.  Significant  features  have  been  labeled  on  the 
plot. 


nence  of  flexural  waves  in  the  side  excitation  of  the  beam,  are  obvious.  However,  a 
substantial  difference  can  be  seen  in  the  testing  in  testing  of  different  excitation  types 
and  mountings. 

In  the  polycarbonate  samples,  the  hung  mounting  showed  fewer  wavetypes  than 
either  of  the  foam  mounting  conditions,  due  to  the  lack  of  restriction  on  the  normal 
vibration  in  the  system.  Under  higher  impedance  mountings,  the  end  excitation 
showed  significantly  more  wavetypes  than  the  side  excitation,  even  some  that  ought 
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to  be  inherent  to  the  side  excitation  (such  as  torsional  waves).  While  waves  relating 
to  the  flexural  motion  of  the  beam  and  longitudinal  vibration  can  easily  be  discerned, 
there  does  not  appear  to  be  any  features  corresponding  to  shear  waves  through  the 
beam. 

The  surrogate  material  samples  showed  similar  characteristics  to  the  polycarbon¬ 
ate  samples;  flexural  waves  and  longitudinal  waves  can  be  seen  in  each  sample,  al¬ 
though  the  75%  flexural  waves  may  be  too  close  in  speed  to  separate  into  normal  and 
transverse  components.  The  50%  sample  showed  significant  damping,  evident  as  a 
one-sided  wavenumber-frequency  plot  that  shows  little  energy  propagation  towards 
the  excitation  location.  The  addition  of  a  stiff  polycarbonate  backing  made  a  sub¬ 
stantial  difference  in  the  50%  sample  case,  masking  all  relevant  wavetypes  that  can 
clearly  be  seen  in  the  foam-mounted  condition.  These  effects  were  not  seen  in  the 
stiffer  75%  sample,  which  is  of  equivalent  density  but  much  stiffen 

The  testing  here  will  serve  as  a  suitable  assessment  of  wave  propagation  properties 
in  the  unknown  surrogate  material;  of  particular  note  is  that  the  50%  material  has 
a  longitudinal  wave  speed  slower  than  that  in  air.  The  subsonic  region  of  wave 
propagation  in  air  will  therefore  be  supersonic  in  the  surrogate,  minimizing  the  effects 
of  evanescent  waves  in  air  in  its  excitation. 
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CHAPTER  7.  CONCLUSIONS 

Evanescent  waves  have  been  shown  to  improve  sound  transmission  across  the  fluid- 
solid  interface  in  subsonic  regions  of  sound  transmission.  By  using  evanescent  waves 
to  reduce  the  refraction  effects  associated  with  large  impedance  difference  boundaries, 
energy  can  be  transmitted  into  solid  materials  via  acoustical  excitation  that  could  not 
be  propagated  in  homogeneous  energy  transmission. 

Evanescent  waves  are  modeled  for  this  work  by  using  a  complex  angle  of  propaga¬ 
tion,  which  produces  a  wavetype  that  decays  normal  to  the  direction  of  propagation; 
this  wavetype  obeys  the  wave  equation,  and  the  complex  angle  can  be  easily  inte¬ 
grated  into  existing  models  of  fluid-solid  interaction.  While  such  waves  would  be 
difficult  to  generate  in  open  space,  the  pressure  and  velocity  of  such  a  wave  can  be 
generated  over  a  limited  space  along  a  surface.  The  complex  angle  of  propagation, 
then,  is  a  convenient  way  to  model  evanescent  waves  in  classical  wave-propagation 
equations. 

The  use  of  evanescent  waves  has  no  beneficial  effect  on  pressure  transmission 
levels,  and  will  make  no  contribution  to  transmission  in  supersonic  regions  of  fluid- 
solid  interaction.  However,  the  intensity  propagation  into  a  material,  particularly 
in  the  subsonic  region,  is  improved  by  the  use  of  incident  evanescent  waves.  At 
incident  angles  above  the  critical  angle,  where  the  trace  wave  speed  on  the  surface 
falls  below  that  of  the  material’s,  refraction  typically  causes  no  energy  propagation 
into  the  material.  The  refraction  of  an  evanescent  wave  does  not  approach  the  critical 
angle,  instead  increasing  the  decay  rate  of  the  wave  into  the  material.  However,  in 
appropriate  situations,  such  a  decay  is  worthwhile  to  allow  for  propagation  into  the 
material. 

Monopoles  have  significant  evanescent  components,  which  can  allow  for  exper¬ 
imental  verification  of  evanescent  wave  effects.  The  decay  of  a  monopole’s  trace 
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pressure  distribution  on  a  surface  decays  with  distance,  allowing  for  a  controlled  ratio 
of  evanescent  to  homogeneous  wave  components.  Using  a  monopole  source  to  judge 
the  effects  of  evanescent  waves  on  fluid-solid  energy  transmission  showed  significant 
contributions  of  the  evanescent  waves  to  the  total  sound  power  transmitted  into  the 
solid. 

While  the  sound  field  of  an  evanescent  wave  would  be  difficult  to  generate,  ap¬ 
proximating  the  boundary  conditions  seen  by  an  evanescent  wave  on  the  fluid-solid 
interaction  surface  would  allow  for  proper  transmission  into  the  material.  A  key  effect 
of  the  evanescent  wave  is  the  phase  difference  between  pressure  and  velocity;  while 
classical  plane  waves  have  in-phase  pressure  and  velocity,  the  addition  of  the  complex 
angle  component  increases  the  phase  difference  towards  a  90  degree  difference. 

Approximating  the  pressure  distributions  made  by  evanescent  waves  on  the  surface 
of  a  material  can  be  done  using  multiple  simple  sources.  A  formulation  is  developed 
that  computes  the  amplitude  and  phase  of  a  monopole  or  dipole  array  based  on 
a  least-squares  approximation  of  an  ideal  pressure  distribution  on  a  surface.  An 
exploration  of  the  parameters  of  the  model  shows  a  heavy  emphasis  on  the  incident 
angle  and  frequency  of  the  incident  wave  being  modeled,  with  angle  of  incidence 
being  a  particularly  important  parameter.  The  geometry  of  the  source  array  shows  a 
proportional  relationship  between  the  standoff  of  the  array  and  the  required  length  to 
obtain  a  good  fit,  with  increased  trace  wavelength  requiring  additional  source  length. 

Attempting  to  approximate  evanescent  waves  on  the  surface  showed  a  significant 
reduction  in  the  quality  of  the  surface  velocity  fit,  most  likely  due  to  the  phase 
difference  between  the  pressure  and  velocity  that  cannot  be  easily  duplicated  using 
simple  sources.  Both  monopole  and  dipole  sources  were  found  to  model  the  desired 
incident  waves  equally  well;  however,  the  dipole  sources  will  have  less  efficiency  due 
to  the  reduced  radiation  efficiency  of  the  dipole. 

Modeling  wave  propagation  through  multi-layer  materials  allows  for  visualization 
of  the  effects  of  geometry  and  material  coupling  on  the  overall  vibration  characteristics 
of  the  system.  A  model  for  propagation  was  formulated  using  the  wave  potentials  in 
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each  layer,  equating  the  material  states  at  each  interface  to  calculate  the  amplitudes 
of  each  potential.  The  use  of  wave  potentials  allows  for  inter-layer  state  calculation, 
as  well  as  intensity  calculation  to  allow  for  measurement  of  energy  propagation  into 
the  material. 

Simulations  from  the  wave  potential  model  show  significant  effects  of  flexural  and 
dilatational  waves,  which  are  formed  from  coupling  of  the  longitudinal  and  shear 
waves.  The  addition  of  thin,  high-impedance  materials  as  barrier  materials  will  not 
affect  these  waves,  as  their  increased  material  stiffness  is  overcome  by  the  low  flexural 
stiffness  resulting  from  their  narrow  width.  However,  removing  the  direct  coupling 
between  the  barrier  and  energetic  layer  causes  an  increase  in  energy  transmission  loss 
into  the  material.  The  effect  of  evanescent  waves  are  put  into  context  by  the  multi¬ 
layer  model;  while  they  can  increase  energy  transmission  in  the  subsonic  region,  other 
avenues  of  energy  propagation  may  provide  higher  energy  transmission,  particularly 
in  lower-speed  materials. 

Measuring  wavetypes  experimentally  is  a  process  dependent  on  careful  planning 
in  order  to  properly  visualize  the  wavetypes  in  the  system.  The  effects  of  mounting 
conditions  and  excitation  location  are  shown  through  the  classification  of  wave  prop¬ 
agation  through  a  polycarbonate  beam  and  two  beams  of  surrogate  materials  with 
unknown  material  properties.  The  polycarbonate  beam’s  properties  are  well-known, 
and  using  the  resulting  vibration  will  allow  a  basis  to  analyze  the  properties  of  the 
unknown  surrogate  materials.  For  the  mounting  conditions,  the  effects  of  any  restric¬ 
tion  to  the  vibration  caused  by  mounting  can  be  seen  in  the  relationship  between  the 
two  flexural  waves  in  the  material;  while  the  normal  motion  dominates  in  the  low- 
impedance  hung  condition,  the  effects  of  the  transverse  wave  on  normal  motion  can 
be  clearly  under  additional  constraint.  When  lower-impedance  materials  are  mounted 
to  materials  with  significantly  higher  impedances,  the  wavetypes  inherent  to  the  ma¬ 
terial  are  masked  by  the  effects  of  the  mounting  materials.  The  excitation  location 
plays  a  significant  role  on  waves  excited  in  the  material;  while  the  side  excitation 


122 


provides  energy  into  the  normal  flexural  wavetype,  exciting  the  materials  from  the 
end  provides  excitation  of  significantly  more  wavetypes. 

The  research  presented  here  provides  an  initial  foray  into  the  theory  and  practical 
considerations  of  evanescent  wave  propagation  into  multi-layer  materials.  Continuing 
upon  this  research  will  allow  for  more  precise  applications  of  the  theory,  especially  in 
regards  to  development  of  novel  excitation  and  detection  technologies. 

7.1  Future  Work 

Of  primary  concern  to  the  research  undertaken  is  the  generation  of  evanescent 
waves;  while  the  theory  presented  here  is  promising,  an  experimental  method  to 
generate  such  waves  will  be  vital  for  exciting  the  systems  of  interest.  The  generation  of 
controlled  pressure  distributions  would,  on  their  own,  be  a  significant  area  of  research, 
it  would  also  allow  for  experimental  measurements  of  evanescent  wave  propagation 
across  fluid-solid  boundaries. 

The  model  developed  here  allows  for  angular  variation  of  source  output,  as  shown 
with  the  dipole  model.  Modifying  the  model  slightly  would  allow  for  arbitrary  angular 
variation  of  each  source,  and  characteristics  representing  real  speakers  could  be  input 
and  calculated.  Additional  modifications  to  the  geometry  can  also  be  made  from 
the  limited  one  explored  here;  for  example;  changing  the  angle  of  the  linear  source 
array,  or  using  creative  geometries.  Optimization  of  geometry  for  a  desired  pressure 
distribution  would  be  a  more  complex  task,  possibly  requiring  the  use  of  a  genetic 
algorithm.  Additional  optimization  of  the  model  may  also  be  possible  by  varying 
the  angle  of  each  speaker  in  the  array.  The  final  goal  of  this  modeling  would  be 
the  creation  and  experimental  verification  of  an  array.  While  other  researchers  have 
developed  arrays  that  allow  for  phase  variation  of  the  sources,  the  use  of  a  64-channel 
output  would  be  the  basis  of  a  setup  that  varies  both  amplitude  and  phase  for  a  given 
frequency. 
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Understanding  of  the  dynamic  effects  of  the  pressure-velocity  phase  difference, 
which  is  an  inherent  property  in  evanescent  waves,  may  provide  additional  avenues 
of  energy  transmission.  In  the  absence  of  evanescent  wave  generation,  additional 
methods  to  excite  evanescent  wave  motion  through  vibrating  a  structure  may  also  be 
explored,  such  as  monopole-likc  radiation  in  a  solid  resulting  from  a  point  excitation. 

Once  evanescent  waves  can  be  found  to  be  generated  in  a  solid,  the  effects  of 
wavetypes  with  concentrated  energy  in  a  particular  frequency  should  be  explored, 
particularly  in  multi-layer  materials.  For  example,  the  ability  to  excite  resonance  in 
a  low-impedance  interior  layer  bounded  by  high-impedance  materials  would  show  the 
transmission  effects  of  evanescent  waves  and  their  ability  to  probe  materials. 

A  detailed  accounting  for  test  methodology  of  unknown  materials  should  be  de¬ 
veloped,  taking  into  account  the  excitation  type  and  mounting  condition.  Several 
wavetypes  of  interest,  such  as  the  shear  wave  and  dilatational  waves,  could  not  be 
seen  or  verified  during  our  tests;  a  novel  excitation  or  measurement  method  may  allow 
for  their  detection  and  classification. 
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APPENDIX  A:  MULTI-LAYER  PROPAGATION  CODE 

The  code  below  comprises  the  multi-layer  propagation  model  developed  in  Chapter 
5.  An  input  structure,  comprising  the  geometry  and  material  properties  of  the  system, 
are  input  along  with  the  frequency,  incident  angle,  and  requested  positions  for  state 
calculation.  The  resulting  output  is  the  wave  potential  amplitudes,  the  states,  and 
resulting  matrices  used  in  calculation. 

The  bulk  of  the  code  comprises  formation  of  the  multi-layer  matrix,  accounting 
for  different  interface  types  and  material  states,  along  with  separate  code  for  the  first 
and  last  layers  (which  contain  distinct  wave  features).  This  code  calls  on  the  code  in 
Appendix  B  for  calculation  of  postion-based  states. 

function  [PSI, states, A, B]  =  MLpot (input , f , inc_angle ,x_req,makeplot) 


if  nargin<5 

raakeplot  =  Jn’ ; 

end 

y„computes  potential  amplitudes  (vector  PSI)  for  multi-layer  fluid-solid 
°/„system  by  constructing  a  linear  system  of  equations  governing  the 
“/interactions  of  each  layer 

“/OUTPUTS : 

“/„  PSI  -  structure  with  of  potential  amplitudes  and  angles 

“/  contains  1  incident  reflected  amplitude,  1  transmitted  output, 

“/  and  4  states  for  each  interior  layer 
“/„  states  -  structure  of  position-dependent  material  states 
“/„  A  -  constructed  matrix 
“/„  B  -  constructed  matrix  in  A*PSI=B 
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% INPUTS : 

1  input  -  structure  describing  materials  and  geometry  of  system 
%  input  structure  (each  layer  is  an  input  in  a  row  vector) : 

°/„  type  -  ’  f’  for  fluid  and  ’s’  for  solid 

°/„  c  -  longitudinal  wavespeed 

°/„  b  -  shear  wavespeed  (’O'  in  fluids) 

%  rho  -  density 

°/„  L  -  length  (’Inf’  on  both  sides) 

%  f  -  frequency 

°/„  inc_angle  -  complex  incident  angle  (degrees) 

%  x_req  -  position  vector  for  state  calculation 
°/„  makeplot  -  make  plot  of  position-dependent  states 

°/„convert  incident  angle  to  radians 

inc_theta  =  real (inc_angle) *pi/180+li*imag(inc_angle) ; 
omega  =  2*pi  *  f ;  °/0f  requency  in  radians 

PSI.f  =  f; 

N  =  length  (input  .type)  ;  °/„number  of  layers 

PSI_N  =  4*(N-2)+2;  °/„number  of  points  in  PSI 

/(calculate  lambda  and  mu  for  each  layer 

mu  =  input . b . “ 2 . * input . rho ; 

lambda  =  input . c . ~2 . * input . rho-2*mu; 


kt  =  omega/input .  c  (1)  *sin(inc_theta)  ;  /(trace  wavelength  -  constant! 


%  material  wavenumbers 
k  =  omega. /input . c ; 
kappa  =  omega . / input . b ; 

°/„calculate  normal  wavenumbers  and  angles  in  each  layer 
k_z  =  sqrt (k . ~2-kt~2) ; 

kappa_z  =  sqrt (kappa. ~2-kt~2) ;  %check  for  imaginary  numbers 
theta  =  asin(kt./k); 
gamma  =  asin(kt . /kappa) ; 

PSI. theta  =  theta; 

PS I. gamma  =  gamma; 

°/„calculate  z-coordinates  of  each  interface 

z=zeros(l ,N-1) ; 

z(l)  =  0; 

for  i  =  2 : N-l ; 

z(i)  =  sum(input . L(2 : i) ) ; 

end 

PSI.z  =  z; 

%  PSI  format: 

%  (1)  incoming  long 
%  (2)  incoming  shear 
%  (3)  outgoing  long 
°/„  (4)  outgoing  shear 


A  =  zeros(PSI_N,PSI_N) ; 
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B  =  zeros(PSI_N, 1) ; 

y„construct  matrix 

°/0first  layer  -  incident  fluid  only,  only  part  with  positive  B  input 
if  input  .type  (2)  ==  ’f’  °/„only  pressure,  velocity  conditions 
"/pressure  condition 

A(l,l)  =  -li*omega* input . rho (1) *exp(li*k_z (1) *-z(l) ) ;  "/reflected  long 
A(l,2)  =  -1* (-li*omega*input . rho (2) . . . 

*exp(-li*k_z(2)*-z(l))) ;  "/transmitted  long 
"/A (1,4)  components  in  "if"  loop  below 

B ( 1 )  =  -1* (-li*omega* input . rho (1) *exp(-li*k_z(l) *-z(l) ) )  ;  "/INCOMING  WAVE 

"/velocity  condition 

A(2,l)  =  li*k(l)*cos(theta(l))*exp(li*k_z(l)*-z(l)) ;  "/reflected  long 
A(2,2)  =  -l*(-li*k(2)*cos(theta(2) )* . . . 

exp(-li*k_z(2)*-z(l))) ;  "/transmitted  long 
"/A (2, 4)  component  in  "if"  loop  below 
B(2)  =  -l*(-li*k(l)*cos(theta(l) )* . . . 

exp(-li*k_z(l)*-z(l))) ;  "/INCOMING  WAVE 

if  N>2  "/incorporate  back-going  waves 

A(l,4)  =  -1* (-li*omega* input . rho (2) * . . . 

exp(li*k_z(2)*-z(l)))  ;  "/back-going  long  (in  second  fluid)  "/for  N>2 

A(2,4)  =  -1* (li*k(2) *cos (theta(2) ) * .  .  . 

exp(li*k_z(2)*-z(l)))  ;  "/back-going  long  for  N>2 


°/„set  shear  coefficients  of  layer  2  to  zero 
A(3,3)  =  1; 
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A(4,5)  =  1; 

end 

r  =  5; 

elseif  input .type (2)  ==’s) 

%stress  relationship 

A(l,l)  =  -1*(  -li  *  omega  *  input. rho(l)  *... 

exp(li  *  -z(l)  *  k_z(l)  ));  °/0reflected  long 
A(1 , 2)  =  (  lambda (2)  *  -k(2)~2  +  2  *  mu(2)  *  -k(2)"2  *... 

cos(theta(2))~2)/(-li*omega)*exp(-li*-z(l)*k_z(2)  );  %transmitted  long 
A(l,3)  =  (2*mu(2) *kappa(2) ~2*cos (gamma(2) ) *sin(gamma(2) ) ) / . . . 
(-li*omega)*exp(-li*-z(l)*kappa_z(2)  );  °/0transmitted  shear 

B ( 1 )  =  (-li*omega*input . rho (1) ) *exp(-li*z(l) *k_z(l) ) ;  %INC0MING  WAVE 

%shear  relationship 
A(2,l)  =  0;  °/„reflected  long 

A(2,2)  =  mu(2) * (2*k(2) ~2*cos (theta(2) ) *sin(theta(2) ) ) / . . . 

(-li*omega)*exp(-li*-z(l)*k_z(2) )  ;  °/„transmitted  long 
A(2,3)  =  mu(2) * (-kappa(2) ~2* (sin(gamma(2) ) ~2-cos (gamma(2) ) ~2) ) / . . . 
(-li*omega)*exp(-li*-z(l)*kappa_z(2)) ;  "/transmitted  shear 

B (2)  =0;  "/INCOMING  WAVE 

"/normal  velocity 

A(3, 1)  =  -1*(  li*k(l)*cos (theta (1) ) *exp(  li*-z(l)*k_z(l))) ; 

A(3,2)  =  -li*k(2)*cos (theta (2) ) *exp(-li*-z (1) *k_z(2) ) ; 
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A(3,3)  =  li*kappa(2)*sin(gamma(2))*exp(-li*-z(l)*kappa_z(2)) ; 

B(3)  =  -li*k(l)*cos(theta(l))*exp(-li*-z(l)*k_z(l)) ; 

if  N>2  °/0add  backward-going  components 

A(l,4)  =  (lambda(2) *  -k(2)~2  +  2*mu(2) *-k(2) ~2* . . . 
cos (theta (2) ) ~2)/ (-li*omega) * . . . 
exp(li*-z(l)*k_z(2)) ; °/0back-going  long 
A(l,5)  =  (2*mu(2) *-kappa(2) ~2*cos (gamma(2) ) *sin(gamma(2) ) )/ . . . 
(-li*omega)  *exp ( li*-z ( 1)  *kappa_z (2) ) ;  °/0back-going  shear 

A(2,4)  =  mu(2) * (2*-k(2) ~2*cos (theta(2) ) *sin(theta(2) ) )/ . . . 

(-li*omega)*exp(li*-z(l)*k_z(2))  ;°/0back-going  long 
A(2,5)  =  mu(2) * (-kappa(2) ~2* (sin(gamma(2) ) ~2-cos (gamma(2) ) ~2) )/ . . . 
(-li*omega)  *exp  ( li*-z  ( 1) *kappa_z  (2) ) ;  °/0back-going  shear 

A(3,4)  =  li*k(2)*cos(theta(2))*. . . 

exp(li*z(l)*k_z(2))  ; °/„back-going  long 
A(3,5)  =  li*kappa(2) *sin(gamma(2) ) * . . . 

exp(li*-z(l) *kappa_z(2) )  ;°/0back-going  shear 

elseif  N==2  7„zero  out  backward  components 
A(4,4)  =  1; 

A(5,5)  =  1; 


B (4)  =  0; 
B (5)  =  0; 
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end 

r  =  4; 

end 

if  N>3 

for  i  =  3:N-1 

if  input .type (i)  ==  ’ f ’ 

if  input .type ( i— 1 ) == J f 5 
°/„pressure 

A(r, l+4*(i-3)+l)  =  -li*omega*input . rho ( i— 1 )  *... 

exp(-  li  *  k_z(i-l)  *  -z(i-l));  °/0incoming  long 
A(r, l+4*(i-3)+3)  =  -li*omega*input . rho (i-1)  *... 

exp(  li  *  k_z(i-l)  *  -z(i-l));  °/„reflected  long 
A(r, l+4*(i-2)+l)  =  -1* (-li*omega*input . rho(i)  *  ... 

exp(-  li  *  k_z(i)  *  -z(i-l)));  "/transmitted  long 
A(r, l+4*(i-2)+3)  =  -1* (-li*omega*input . rho (i)  *  ... 
exp(  li  *  k_z(i)  *  — z (i— 1) ) ) ;  "/back-going  long 

"/velocity 

A(r+1 , 1+4* ( i— 3) +1 )  =  -k(i-l) *cos (theta(i-l) ) * . . . 

exp(-li  *  k_z(i-l)  *  -z(i-l)); 

A(r+1 , 1+4* (i-3)+3)  =  k(i-l) *cos (theta(i-l) ) * . . . 

exp(li  *  k_z(i-l)  *  -z(i-l)); 

A(r+1 , 1+4* (i-2)+l)  =  -1* (-k(i) *cos (theta(i) ) * . . . 

exp(-li  *  k_z(i)  *  -z(i-l))); 

A(r+1 , 1+4* (i-2)+3)  =  -1* (k(i) *cos (theta(i) ) * . . . 

exp(li  *  k_z(i)  *  -z(i-l))); 


°/„shear  components  in  i-th  fluid  =  0 

A(r+2 , 1+4* (i-2)+2)  =  1; 

A(r+3,l+4*(i-2)+4)  =  1; 

r  =  r+4; 

elseif  input .type ( i— 1 ) ==  ^ s ’  %solid-f luid  interface 

°/„normal  stress 

A(r, l+4*(i-3)+l)  =  (lambda(i-l) *-k(i-l) ~2+ . . . 
2*mu(i-l)*-k(i-l)~2*cos(theta(i-l))~2)*. . . 
exp(-li*k_z(i-l)*-z(i-l) )/(-li*omega) ; 

A(r, l+4*(i-3)+2)  =  2*mu(i-l) *kappa(i-l) ~2* . .  . 
cos (gamma(i-l) ) *sin(gamma(i-l) ) * . . . 
exp(-li*kappa_z(i-l)*-z(i-l))/ (-li*omega) ; 

A(r, l+4*(i-3)+3)  =  (lambda(i-l) *-k(i-l) ~2+ . .  . 
2*mu(i-l) *-k(i-l) ~2*cos (theta(i-l) ) "2) * . . . 
exp(li*k_z(i-l)*-z(i-l))/ (-li*omega) ; 

A(r, l+4*(i-3)+4)  =  2*mu(i-l) *-kappa(i-l) ~2* .  .  . 
cos (gamma(i-l) ) *sin(gamma(i-l) ) * . . . 
exp(li*kappa_z(i-l)*-z(i-l))/ (-li*omega) ; 

A(r, l+4*(i-2)+l)  =  -1* (-li*omega*input . rho(i) * 
exp(-li*k_z(i)*-z(i-l))) ; 

A(r, l+4*(i-2)+3)  =  -1* (-li*omega*input . rho (i) * 
exp(li*k_z(i)*-z(i-l) ) ) ; 

°/„transverse  stress 

A(r+1 , 1+4* (i-3)+l)  =  mu(i-l)*(2*k(i-l)~2*. . . 


cos (theta(i-l) ) *sin(theta(i-l) ) ) * . . . 

exp(-li*k_z(i-l)*-z(i-l) )/(-li*omega) ; 

A(r+1 , 1+4* (i-3)+2)  =  mu(i-l) * (-kappa(i-l) "2* . . . 
(sin(gamma(i-l) ) ~2-cos (gamma (i-1) ) “2) ) * . . . 
(exp(-li*kappa_z(i-l)*-z(i-l)))/ (-li*omega) ; 

A(r+1 , 1+4* (i-3)+3)  =  mu(i-l) *(2*-k(i-l) ~2* . . . 
cos (theta(i-l) ) *sin(theta(i-l) ) ) * . . . 
(exp(li*k_z(i-l)*-z(i-l)))/(-li*omega) ; 

A(r+1 , 1+4* (i-3)+4)  =  mu(i-l) * (-kappa(i-l) "2* . . . 
(sin(gamma(i-l) ) "2-cos (gamma (i-1) ) ~2) )* . . . 
exp(li*kappa_z(i-l)*-z(i-l))/ (-li*omega) ; 

A(r+1 , 1+4* (i-2)+l)  =  0; 

A(r+1 , 1+4* (i-2)+3)  =  0; 

7„norraal  velocity 

A(r+2 , 1+4* ( i— 3) +1 )  =  -li*k(i-l) *cos (theta(i-l) ) * . . 
exp(-li*k_z(i-l)*-z(i-l)) ; 

A(r+2 , 1+4* (i-3)+2)  =  li*kappa(i-l) *sin(gamma(i-l) ) 
exp(-li*kappa_z(i-l)*-z(i-l))  ; 

A(r+2 , 1+4* (i-3)+3)  =  li*k(i-l) *cos (theta(i-l) ) * . . . 
exp(li*k_z(i-l)*-z(i-l)) ; 

A(r+2, l+4*(i-3)+4)  =  li*kappa(i-l) *sin(gamraa(i-l) ) 
exp(li*kappa_z(i-l)*-z(i-l)) ; 

A(r+2 , 1+4* ( i— 2) +1 )  =  -1* (— li*k(i) *cos (theta(i) ) * . . 
exp(-li*k_z(i)*-z(i-l)))  ; 

A(r+2 , 1+4* (i-2)+3)  =  -1* (li*k(i) *cos (theta(i) ) * . . . 
exp(li*k_z(i)*-z(i-l))) ; 


*/.  set  shear  components  in  i-th  fluid  to  0 


A(r+3, l+4*(i-2)+2)  =  1; 
A(r+4, l+4*(i-2)+4)  =  1; 

r  =  r+5; 

end 

elseif  input .  typed)  ==  J  s’ 
if  input .type ( i— 1 ) ==  ^ s’ 


°/„normal  stress 

A(r, l+4*(i-3)+l)  =  (lambda(i-l) *-k(i-l) ~2+ . . . 
2*mu(i-l) *-k(i-l) ~2*cos (theta(i-l) ) ~2)  * . . 
exp(-li*k_z(i-l)*-z(i-l) )/(-li*omega) ; 
A(r, l+4*(i-3)+2)  =  2*mu(i-l) *kappa(i-l) ~2* . . . 
cos (gamma (i-1) ) *sin(gamma(i-l) ) * . . . 
exp(-li*kappa_z(i-l) *-z(i-l) ) / (-li*omega) 
A(r, l+4*(i-3)+3)  =  (lambda(i-l) *-k(i-l) “2+ . . . 
2*mu(i-l) *-k(i-l) ~2*cos (theta (i-1) ) ~2)  * . . 
exp(li*k_z(i-l)*-z(i-l))/ (-li*omega) ; 

A(r, l+4*(i-3)+4)  =  2*mu(i-l) *-kappa(i-l) ~2* . . 
cos (gamma (i-1) ) *sin(gamma(i-l) ) * . . . 
exp(li*kappa_z(i-l)*-z(i-l))/ (-li*omega) ; 
A(r, l+4*(i-2)+l)  =  -1* ( (lambda(i) *-k(i) ~2+ . . . 
2*mu(i) *-k(i) ~2*cos (theta (i) ) "2) * . . . 
exp(-li*k_z(i)*-z(i-l) )/(-li*omega) ) ; 

A(r, l+4*(i-2)+2)  =  -1* (2*mu(i) *kappa(i) ~2* . . . 
cos (gamma(i) ) *sin(gamma(i) ) * . . . 
exp(-li*kappa_z(i)*-z(i-l))/ (-li*omega) ) ; 
A(r, l+4*(i-2)+3)  =  -1* ( (larabda(i) *-k(i) "2+ . . . 
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2*mu(i) *-k(i) ~2*cos (theta(i) ) ~2) * . . . 
exp(li*k_z(i)*-z(i-l))/ (-li*omega) ) ; 

A(r,l+4*(i-2)+4)  =  -1* (2*mu(i) *-kappa(i) ~2* . . . 
cos  (gamma(i) )  *sin(gamina(i) )  * . . . 
exp(li*kappa_z(i)*-z(i-l))/ (— li Iomega)) ; 

°/„tangential  stress 

A(r+l,l+4*(i-3)+l)  =  mu(i-l)*(2*k(i-l)~2*. . . 
cos (theta(i-l) ) *sin(theta(i-l) ) ) * . . . 
exp(-li*k_z(i-l)*-z(i-l))/ (-li*omega) ; 

A(r+1 , 1+4* (i-3)+2)  =  mu(i-l) * (kappa (i-1) ~2* . . . 
(-sin(gamma(i-l) ) ~2+cos (gamma(i-l) ) ~2) ) * . . . 
(exp(-li*kappa_z(i-l)*-z(i-l)))/ (-li*omega) ; 

A(r+1 , 1+4* (i-3)+3)  =  mu(i-l)*(2*-k(i-l)~2*. . . 
cos (theta(i-l) ) *sin(theta(i-l) ) ) * . . . 

(exp(li*k_z(i-l)*-z(i-l)))/(-li*omega) 5 

A(r+1 , 1+4* (i-3)+4)  =  mu(i-l) * (kappa (i-1) ~2* . . . 
(-sin(gamma(i-l) ) ~2+cos (gamma(i-l) ) ~2) ) * . . . 
exp(li*kappa_z(i-l)*-z(i-l))/ (-li*omega) ; 

A(r+1 , 1+4* (i-2)+l)  =  -l*(mu(i)*(2*k(i)~2*. . . 
cos (theta(i) ) *sin(theta(i) ) ) * . . . 
exp(-li*k_z(i)*-z(i-l) )/(-li*omega) ) ; 

A(r+1 , 1+4* (i-2)+2)  =  -l*(mu(i) *(kappa(i) “2* . . . 
(-sin(gamma(i) ) ~2+cos (gamma(i) ) ~2) ) * . . . 
(exp(-li*kappa_z(i)*-z(i-l)))/ (-li*omega) ) ; 

A(r+1 , 1+4* (i-2)+3)  =  -1* (mu(i) * (2*-k(i) ~2* . . . 
cos (theta(i) ) *sin(theta(i) ) ) * . . . 
(exp(li*k_z(i)*-z(i-l)))/ (-li*omega) ) ; 

A(r+1 , 1+4* (i-2)+4)  =  -1* (mu(i) * (kappa(i) ~2* . . . 


(-sin(gamma(i) ) ~2+cos (gamma(i) ) ~2) ) * . . . 
exp(li*kappa_z(i)*-z(i-l))/ (-li*omega)) ; 


141 


y„normal  velocity 

A(r+2 , 1+4* (i-3)+l)  =  li*  -k(i-l) *cos (theta(i-l) ) * . . . 
exp(-li*k_z(i-l)*-z(i-l)) ; 

A(r+2 , 1+4* (i-3)+2)  =  li*kappa(i-l) *sin(gamraa(i-l) ) * . . . 
exp(-li*kappa_z(i-l)*-z(i-l)) ; 

A(r+2 , 1+4* (i-3)+3)  =  li*k(i-l) *cos (theta(i-l) ) * . . . 
exp(li*k_z(i-l)*-z(i-l)) ; 

A(r+2, l+4*(i-3)+4)  =  li*kappa(i-l)*sin(gamma(i-l) )* . . . 
exp(li*kappa_z(i-l)*-z(i-l)) ; 

A(r+2 , 1+4* ( i— 2) +1 )  =  -1* ( 1 i * — k ( i ) *cos (theta(i) ) * . . . 
exp(-li*k_z(i)*-z(i-l))) ; 

A(r+2 , 1+4* (i-2)+2)  =  -1* (li*kappa(i) * s in (gamma (i) ) * . . . 
exp(-li*kappa_z(i)*-z(i-l))) ; 

A(r+2 , 1+4* (i-2)+3)  =  -1* (li*k(i) *cos (theta(i) ) * . . . 
exp(li*k_z(i)*-z(i-l))) ; 

A(r+2, l+4*(i-2)+4)  =  -1* (li*kappa(i) *sin(gamma(i) ) * . . . 
exp(li*kappa_z(i)*-z(i-l))) ; 

%tangential  velocity 

A(r+3, l+4*(i-3)+l)  =  li*k(i-l) *sin(theta(i-l) ) * . . . 
exp(li  *  -k_z(i-l)  *  -z(i-l)); 

A(r+3, l+4*(i-3)+2)  =  -li*-kappa(i-l) *cos (gamma(i-l) ) * . 
exp(li  *  -kappa_z(i-l)  *  -z(i-l)); 

A(r+3, l+4*(i-3)+3)  =  li*k(i-l) *sin(theta(i-l) ) * . . . 
exp(li  *  k_z(i-l)  *  -z(i-l)); 

A(r+3, l+4*(i-3)+4)  =  -li*kappa(i-l) *cos (gamma(i-l) ) * . . 
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exp(li  *  kappa_z(i-l)  *  -z(i-l)); 

A(r+3, l+4*(i-2)+l)  =  -1*(  li*k(i) *sin(thetad) ) * . . . 

exp(li  *  -k_z(i)  *  -z(i-l)  )  ); 

A(r+3, l+4*(i-2)+2)  =  -1* (-li*-kappa(i) *cos (garama(i) ) * . 

exp(li  *  -kappa_z(i)  *  -z(i-l))  ); 
A(r+3d+4*(i-2)+3)  =  -1*(  li*kd) *sin(thetad) ) * .  .  . 

exp(li  *  k_z(i)  *  -z(i-l))  ); 

A(r+3d+4*(i-2)+4)  =  -1* (-li*kappad) *cos  (gamraa(i) ) *  .  . 
exp(li  *  kappa_z(i)  *  -z(i-l))  ); 


r  =  r+4; 

elseif  input  .type ( i— 1 )  ==  ^ f  ’  °/0f  luid-solid  interface 

7„norraal  stress 

A(rd+4*d-3)+l)  =  -1*  (-li*omega*input .  rho(i-l) * .  .  . 
exp(-li*k_zd-l)*-zd-l) ) ) ; 

A(rd+4*(i-3)+3)  =  -1* (-li*omega*input .  rho(i-l) * .  .  . 
expdi*k_zd-l)*-zd-l))) ; 

A(rd+4*d-2)+l)  =  (lambda(i) *-kd) ~2+2*mud) * .  .  . 
-kd)  ~2*cos  (thetad) )  ~2)  * . . . 
exp(-li*k_zd)*-zd-l)  )/(-li*omega)  ; 

A(r,l+4*(i-2)+2)  =  2*mud) *kappad) ~2* . . . 
cos (gamma(i) ) * s in ( gamma d) ) * . . . 
exp(-li*kappa_zd)*-zd-l))/ (-li*omega) ; 

A(r,l+4*(i-2)+3)  =  (lambda(i) *-kd) ~2+2*mud) * . . . 
-k(i)  ~2*cos  (theta(i) )  “2)  * .  .  . 
expdi*k_zd)*-zd-l))/ (-li*omega) ; 

A(r,l+4*(i_2)+4)  =  2*mud)  *_kappad)  ~2*  .  .  . 
cos (gamma(i) ) *sin(gammad) ) * .  .  . 
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exp(li*kappa_z(i)*-z(i-l))/ (-li*omega) ; 

Ztransverse  stress 

A(r+l,l+4*(i-3)+l)  =  0; 

A(r+1 , 1+4* (i-3)+3)  =  0; 

A(r+1 , 1+4* (i-2)+l)  =  mu(i) *2*  k(i) ~2*cos (theta(i) ) * . . . 
sin (theta (i) ) *exp(-li*k_z(i) *-z(i-l) ) / (—11 *omega) ; 

A(r+1 , 1+4* (i-2)+2)  =  mu(i) * (-kappa(i) "2* . . . 
(sin(gamraa(i) ) "2-cos (gamma(i) ) ~2) )* . . . 
exp(-li*kappa_z(i)*-z(i-l))/ (-li*omega) ; 

A(r+1 , 1+4* (i-2)+3)  =  mu(i) * (2*-k(i) ~2*cos (theta(i) ) * . . 
sin (theta (i) ) ) *exp(li*k_z(i) *-z(i-l) ) / (-li*omega) ; 

A(r+1 , 1+4* (i-2)+4)  =  mu(i) * (-kappa(i) ~2* . . . 
(sin(gamraa(i) ) "2-cos (gamma(i) ) ~2) )* . . . 
exp(li*kappa_z(i)*-z(i-l))/ (-li*omega) ; 

y„normal  velocity 

A(r+2 , 1+4* (i-3)+l)  =  -1* (-li*k(i-l) *cos (theta(i-l) ) * . . 
exp(-li*k_z(i-l)*-z(i-l))) ; 

A(r+2 , 1+4* (i-3)+3)  =  -1* (li*k(i— 1) *cos (theta(i-l) ) * . . . 
exp(li*k_z(i-l)*-z(i-l))) ; 

A(r+2 , 1+4* (i-2)+l)  =  -li*k(i) *cos (theta(i) ) * . . . 
exp(-li*k_z(i)*-z(i-l) ) ; 

A(r+2 , 1+4* ( i— 2) +2)  =  li*kappa(i) *sin (gamma (i) ) * . . . 
exp(-li*kappa_z(i)*-z(i-l)) ; 

A(r+2 , 1+4* (i-2)+3)  =  li*k(i) *cos (theta(i) ) * . . . 
exp(li*k_z(i)*-z(i-l)) ; 

A(r+2, l+4*(i-2)+4)  =  li*kappa(i) *sin (gamma (i) ) * . . . 
exp(li*kappa_z(i)*-z(i-l)) ; 
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r  =  r+3 ; 

end 

end 

end 

end 

if  N>2 

i  =  N;  °/„f inal  layer  -  only  outgoing  long  (into  infinite  fluid) 

if  input .type (N— 1 )  ==  ’f’ 

°/„pressure 

A(PSI_N-l,l+4*(i-3)+l)  =  -li*omega*input.rho(i-l) . . . 

*  exp(-  li  *  k_z(i-l)  *  -z(i-l));  /(incoming  long 
A(PSI_N-1 , l+4*(i-3)+3)  =  -li*omega*input.rho(i-l) . . . 

*  exp(  li  *  k_z(i-l)  *  -z(i-l));  /(reflected  long 
A(PSI_N-1 , 1+4* (i-2)+l)  =  -1* (-li*omega* input . rho (i) . . . . 

*  exp(  -li  *  k_z(i)  *  -z(i-l)));  "/(transmitted  long  long 

°/„velocity 

A(PSI_N,l+4*(i-3)+l)  =  -k(i-l) *cos (theta(i-l) )  .  .  . 

*exp(-li  *  k_z(i-l)  *  -z(i-l)); 

A(PSI_N, l+4*(i-3)+3)  =  k(i-l) *cos (theta(i-l) ) . . . 

*exp(li  *  k_z(i-l)  *  -z(i-l)); 

A(PSI_N, l+4*(i-2)+l)  =  -1* (-k(i) *cos (theta(i) ) . . . 

*exp(-li  *  k_z(i)  *  -z(i-l)));  /(outgoing  long 


elseif  input .type (N— 1)  ==  ’s 
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°/„normal  stress 

A(PSI_N-2,l+4*(i-3)+l)  =  (lambda(i-l) *-k(i-l) ~2+ . . . 

2*mu(i-l) *-k(i-l) ~2*cos (theta(i-l) ) ~2) . . . 

*exp(-li*k_z(i-l) *-z(i-l) )/(-li*omega) ; 

A(PSI_N-2, l+4*(i-3)+2)  =  2*mu(i-l) *kappa(i-l) ~2* . . . 

cos(gamma(i-l))*sin(gamma(i-l)) . . . 

*exp(-li*kappa_z(i-l)*-z(i-l))/ (-li*omega) ; 

A(PSI_N-2,l+4*(i-3)+3)  =  (lambda(i-l) *-k(i-l) “2+ . . . 

2*mu(i-l) *-k(i-l) ~2*cos (theta(i-l) ) "2) .  .  . 
*exp(li*k_z(i-l)*-z(i-l))/ (-li*omega) ; 

A(PSI_N-2, 1+4* (i-3)+4)  =  2*mu(i-l) *-kappa(i-l) ~2* . . . 

cos(gamma(i-l))*sin(gamma(i-l))*exp(li*kappa_z(i-l) . . . 

*-z(i-l) )/ (-li*omega) ; 

A(PSI_N-2,l+4*(i-2)+l)  =  -1* (_li*omega*input . rho (i) * . . . 
exp(-li*k_z(i)*-z(i-l))) ; 

/(transverse  stress 

A(PSI_N-l,l+4*(i-3)+l)  =  mu(i-l)*2*k(i-l) ~2*cos(theta(i-l) ) . . . 
*sin(theta(i-l))*exp(-li*k_z(i-l)*-z(i-l))/ (-li*omega) ; 

A(PSI_N-l,l+4*(i-3)+2)  =  mu(i-l) * (-kappa(i-l) ~2* . . . 

( s in (gamma (i-1) ) ~2-cos (gamma(i-l) ) ~2) ) . . . 

*(exp(-li*kappa_z(i-l)*-z(i-l)))/ (-li*omega) ; 

A(PSI_N-l,l+4*(i-3)+3)  =  mu(i-l)*(2*-k(i-l)~2*. . . 

cos (theta(i-l) ) *sin(theta(i-l) ) ) . . . 

*(exp(li*k_z(i-l)*-z(i-l)))/ (-li*omega) ; 

A(PSI_N-l,l+4*(i-3)+4)  =  mu(i-l) * (-kappa(i-l) ~2* . . . 

(sin(gamma(i-l) ) ~2-cos (gamma(i-l) ) ~2) ) . . . 

*exp(li*kappa_z(i-l)*-z(i-l))/ (-li*omega) ; 


A(PSI_N-l,l+4*(i-2)+l)  =  0; 
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y„normal  velocity 

A(PSI_N,l+4*(i-3)+l)  =  -li*k(i-l) *cos (theta(i-l) ) . . . 
*exp(-li*k_z(i-l)*-z(i-l)) ; 

A(PSI_N, l+4*(i-3)+2)  =  li*kappa(i-l) *sin(gamma(i-l) ) . . . 

*exp(-li*kappa_z(i-l)*-z(i-l)) ; 

A(PSI_N, l+4*(i-3)+3)  =  li*k(i-l) *cos (theta(i-l) ) . . . 
*exp(li*k_z(i-l)*-z(i-l) ) ; 

A(PSI_N, l+4*(i-3)+4)  =  li*kappa(i-l) *sin (gamma (i-1) ) . . . 

*exp(li*kappa_z(i-l)*-z(i-l) ) ; 

A(PSI_N,l+4*(i-2)+l)  =  -l*(-li*k(i)*cos(theta(i) ) . . . 
*exp(-li*k_z(i)*-z(i-l)) ) ; 


end 

end 

PSI .pot  =  inv(A) *B ; 

if  nargin>3 

states  =  MLstate_calc (input,  PSI,  x_req,makeplot) ; 

else 

states  =  ’  ’  ; 


end 
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APPENDIX  B:  MULTI-LAYER  STATE  CALCULATION  CODE 

This  code  calculates  position-based  states  for  the  model  presented  in  Chapter  5. 
When  given  the  input  parameters  for  the  model  along  with  the  solutions  obtained  us¬ 
ing  the  code  in  Appendix  A,  the  calculation  of  pressure,  stress,  velocity,  and  intensity 
is  performed  at  the  requested  positions. 

The  inputs  for  the  model  are  the  input  parameters  of  geometry  and  material  prop¬ 
erties  in  the  system  (the  same  input  used  in  Appendix  A),  the  solutions  obtained  by 
the  code  in  Appendix  A,  and  the  position  vector  for  calculation.  The  resulting  out¬ 
put  is  a  structure  containing  all  the  material  states,  including  intensity  and  reference 
input  parameters,  for  each  position. 

function  states  =  MLst ate (input,  PSI,  x,makeplot) 

“/EDITS  4/30  -  added  conditions  for  infinite  solid  interface 

if  nargin<4 

raakeplot  =  ’n’ ; 

end 

“/compute  states  in  material  (as  a  function  of  incoming  pressure) 

“/initial  propagation  angle  basked  into  solution? 

N  =  length ( input .  type)  ;  “/number  of  layers 

omega  =  2*pi  *  PSI.f; 


mu  =  input  .b .  ~2 .  *  input .  rho ; 
lambda  =  input . c . ~2 . *input . rho-2*mu; 


°/„  compute  x-coordinates  of  borders 
z=zeros(l ,N-1) ; 

z(l)  =  0; 
for  i  =  2 : N-l ; 

z(i)  =  sum(input . L(2 : i) ) ; 

end 

k  =  omega . / input . c ; 

kappa  =  omega . / input . b ; 

k_x  =  k. *sin(PSI .theta) ; 

k_z  =  k. *cos(PSI .theta) ; 

kappa_x  =  kappa. *sin(PSI .gamma) ; 

kappa_z  =  kappa. *cos (PSI .gamma) ; 

•/.if  borders  (z)  not  present,  add  to  system 

for  i  =  l:length(z) 

ind  =  f ind(z (i) >=x, 1 , ’ last ’ ) ; 

%  add  a  second  point  if  already  there 
if  length ( ind) ==0 

x  =  [z(l)  z ( 1 )  x]  ; 

else 

x  =  [x(l:ind)  z(i)  x(ind+l : end)]  ; 
if  x(ind)  ~=  z(i) 

x  =  [x(l:ind)  z(i)  x(ind+l : end)] 
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end 

end 

end 

states. x  =  x; 

"preference  intensity  -  intensity  at  z  =  0; 

I_ref  =  l/2*real (-li*omega*input . rho (1)  *  conj (li*k_z(l)*  (-1))); 

j  =  1;  °/0start  on  first  layer 

for  i  =  1: length (x)  %cycle  through  all  points 

if  j  ==  1 

Zz_in  =  -li*omega*input .rho(l)  *  (l*exp(-li*k_z(l)*-x(i))) ; 
Zz_out  =  -li*omega*input . rho(l)  *  (PSI .pot (1) * . . . 
exp(li*k_z(l)*-x(i))) ; 

Vz_in  =  li*k_z(l)*  (-l*exp(-li*k_z(l) *-x(i) ) ) ; 

Vz_out  =  li*k_z(l)*  (PSI . pot (1) *exp(li*k_z(l) *-x(i) ) ) ; 


states. Zz(i)  =  Zz_in+Zz_out ; 

states. Zx(i)  =  0; 

states. Vz(i)  =  Vz_in+Vz_out ; 

states. Vx(i)  =  li*k_x( j ) * (PSI .pot (1) *exp(li*k_z(l) *-x(i) )+ . . . 
l*exp(-li*k_z(l)*-x(i))) ; 


states . I (1 , i)  =  l/2*real(Zz_out*conj (Vz_out) ) ; 
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elseif  j  ==  N 

if  input .type(j)  ==  5 f’ 

Zz_out  =  -li*omega*input . rho (N)  *... 

(PSI ,pot(l+4*(N-2)+l)*exp(-li*k_z(N)*-x(i))) ; 

Vz_out  =  -li*k_z(N)*PSI.pot(l+4*(N-2)+l)*. . . 
exp(-li*k_z(N)*-x(i) ) ; 

states. Zz(i)  =  Zz_out; 
states. Zx(i)  =  0; 
states. Vz(i)  =  Vz_out; 

states. Vx(i)  =  li*k_x(N) *PSI .pot (1+4* (N-2)+l) * . . . 
exp(-li*k_z (N) *-x(i) ) ; 

states . I (1+4* (N-2)+l , i)  =  l/2*real (Zz_out*conj (Vz_out) ) ; 
elseif  input .type(j)  ==  ’s’ 

Zz_forward  =  ( (lambda( j ) *-k( j ) ~2+2*mu( j ) *-k_z ( j ) ~2) * . . . 

PSI .pot (1+4* ( j-2)+l) *exp(-li*k_z( j ) *-x(i) ) . . . 

+2*mu ( j ) *kappa_z ( j ) *kappa_x ( j ) *PSI . pot ( 1+4* ( j -2) +2) * . . 
exp(-li*kappa_z(j)*-x(i)))/ (-li*omega) ; 

Vz_f orward  =  (li*k_z( j ) *-PSI .pot (1+4* ( j-2)+l) * . . . 
exp(-li*k_z(j)*-x(i))+li*kappa_x(j)*. .  . 

PSI .pot (1+4* ( j  —2) +2) *exp(-li*kappa_z( j ) *-x(i) ) ) ; 

Zx_forward  =  mu( j ) * (2*k_x(j ) *k_z( j ) *PSI .pot (1+4* ( j-2)+l) * . 
exp(-li*k_z ( j ) *-x(i) )+(-kappa_x( j ) ~2+kappa_z( j ) ~2) * . . . 
PSI . pot (1+4* ( j  —2) +2) *exp(-li*kappa_z( j ) *-x(i) ) ) / . . . 


(-li*omega) ; 


Vx_forward  =  (li*k_x( j ) *PSI .pot (1+4* (j-2)+l) * . . . 
exp(-li*k_z ( j ) *-x(i) )-li*kappa_z( j ) * . . . 

-PSI .pot (l+4*(j-2)+2)*exp(-li*kappa_z(j)*-x(i))) 

states. Zz(i)  =  Zz_forward; 
states. Zx(i)  =  Zx_forward; 
states. ¥z(i)  =  Vz_forward; 
states. Vx(i)  =  Vx_forward; 


states . I (1+4* (j-2)+l , i)  =  l/2*real (Zz_f orward* .. . 
conj (Vz_f orward) ) ; 

states . I (1+4* (j-2)+2 , i)  =  l/2*real (Zx_f orward* .. . 
conj (Vx_f orward) ) ; 


end 


else 


if  input .type(j)  ==  ’ 

Zz_forward  =  -li*omega* input . rho ( j ) * . . . 

(PSI.pot(l+4*(j-2)+l)*exp(-li*k_z(j)*-x(i))) ; 
Zz_backward  =  -li*omega*input . rho (j )  *... 

(PSI .pot (1+4* ( j-2)+3) *exp(li*k_z( j )*-x(i) ) ) ; 
Vz_forward  =  li*k_z(j)*. . . 

-PSI .pot (  l+4*(j-2)+l)*exp(-li*k_z(j)*-x(i)) ; 
Vz_backward  =  li*k_z(j)* . . . 


PSI . pot (1+4* ( j -2) +3) *exp ( li*k_z ( j ) *-x (i) ) ; 


states. Zz(i)  =  Zz_forward  +  Zz_backward; 
states. Zx(i)  =  0; 

states. Vz(i)  =  Vz_forward  +  Vz_backward; 
states. Vx(i)  =  li*k_x( j ) * (PSI . pot (  l+4*(j-2)+l)*. . . 
exp(-li*k_z ( j ) *-x(i) )+PSI . pot (1+4* ( j-2)+3) * .  .  . 
exp(li*k_z(j)*-x(i))) ; 

states . I (1+4* ( j  —2) +1 , i)  =  l/2*real (Zz_f orward* . . . 
conj (Vz_f orward) ) ; 

states . I (1+4* (j-2)+3, i)  =  l/2*real (Zz_backward* . . . 
conj (Vz_backward) ) ; 

elseif  input .type(j)  ==  ’s’ 

Zz_forward  =  ( (lambda( j ) *-k( j ) ~2+2*mu( j ) *-k_z ( j ) ~2) * .  . 
PSI .pot (1+4* ( j-2)+l) *exp(-li*k_z( j ) *-x(i) )+2*mu( j ) 
kappa_z(j)*kappa_x(j)*PSI ,pot(l+4*(j-2)+2)* . . . 
exp(-li*kappa_z(j)*-x(i)))/ (-li*omega) ; 

Zz_backward  =( (lambda( j ) *-k( j ) ~2+2*mu( j ) *-k_z ( j ) ~2) * . . 
PSI . pot (1+4* ( j-2)+3) *exp(li*k_z( j ) *-x(i) )+2*mu( j ) * 
-kappa_z ( j ) *kappa_x ( j ) *PSI . pot ( 1+4* ( j -2) +4) * . . . 
exp(li*kappa_z(j)*-x(i)))/ (-li*omega) ; 

Vz_f orward  =  (li*k_z(j)*-PSI ,pot(l+4*(j-2)+l)* . . . 
exp(-li*k_z ( j ) *-x(i) )+li*kappa_x( j ) * . . . 

PSI .pot (1+4* ( j  —2) +2) *exp(-li*kappa_z( j ) *-x(i) ) ) ; 
¥z_backward  =  (li*k_z(j)*PSI ,pot(l+4*(j-2)+3)* . . . 
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exp(li*k_z(j)*-x(i))+li*kappa_x(j)* . . . 

PSI .pot (1+4* ( j-2)+4) *exp(li*kappa_z( j )*-x(i))); 

Zx_forward  =  mu( j ) * (2*k_x( j ) *k_z( j ) *PSI .pot (1+4* ( j-2)+l) * . . . 
exp(-li*k_z(j)*-x(i)) . . . 

+(-kappa_x(j ) ~2+kappa_z( j ) ~2) *PSI . pot (1+4* ( j-2)+2) * . . . 
exp(-li*kappa_z(j)*-x(i)))/ (-li*omega) ; 

Zx_backward  =  mu( j ) * (2*k_x( j ) *k_z(j ) * . . . 

-PSI.pot(l+4*(j-2)+3)*exp(li*k_z(j)*-x(i)) • • • 

+(-kappa_x(j ) ~2+kappa_z( j ) ~2) *PSI . pot (1+4* ( j-2)+4) * . . . 
exp(li*kappa_z(j)*-x(i)))/ (-li*omega) ; 

¥x_forward  =  (li*k_x(j ) *PSI .pot (1+4* (j-2)+l) * .  .  . 
exp(-li*k_z ( j ) *-x(i) )-li*kappa_z( j ) * . . . 

-PSI .pot (1+4* (j-2)+2) *exp(-li*kappa_z(j ) *-x(i) ) ) ; 
¥x_backward  =  (li*k_x(j)*PSI.pot(l+4*(j-2)+3)*. . . 

exp(li*k_z(j)*-x(i))-li*kappa_z(j)* . . . 

PSI .pot (l+4*(j-2)+4)*exp(li*kappa_z(j)*-x(i))) ; 

states. Zz(i)  =  Zz_forward  +  Zz_backward; 

states. Zx(i)  =  Zx_forward  +  Zx_backward; 

states. ¥z(i)  =  ¥z_forward  +  ¥z_backward; 

states. ¥x(i)  =  ¥x_forward  +  ¥x_backward; 

states . I (1+4* (j-2)+l , i)  =  l/2*real (Zz_f orward* . . . 
conj (¥z_f orward) ) ; 

states. I(l+4*(j-2)+2,i)  =  l/2*real (Zx_f orward* .. . 
conj (¥x_f orward) ) ; 

states . I (1+4* (j-2)+3, i)  =  l/2*real (Zz_backward* . . . 


conj (Vz_backward) ) ; 

states . I (1+4* (j-2)+4, i)  =  l/2*real (Zx_backward* . . 
conj (Vx_backward) ) ; 

end 

end 

if  j  ~=N 

if  z(j)  ==  x(i)  &&  z(j)  ==  x(i+l) 

% compute  on  both  sides  of  the  material  border. . 

% (increases  j  by  one  after  first  border  x  value) 

j  =  j+i; 


end 

end 

end 

/(normalize  against  incoming  pressure 

states. Pin  =  - li*omega* input . rho (1) ; 

states. Zz  =  states . Zz/states . Pin; 
states. Zx  =  states . Zx/states . Pin; 
states. Vz  =  states .Vz/states . Pin; 
states. Vx  =  states .Vx/states . Pin; 

states. 1=  states . I . /I_ref ; 
states. I_ref  =  I_ref; 


states. 11(1)  =  states . I (1 , 1)+1 ; 


for  i  =  2:length(x) 

if  states. x(i)<=0  &&  states .x(i-l) ~=0 
states. II  =  states . I (1 )+l ; 

else 

states. I 1 ( i )  =  0; 

end 

end 


“/ocalculat  Iz,  lx  for  each  component 

if  strcmp(makeplot ,  ’n’)~=l 
%  °/„plot  data 
figure 
subplot (221) 
hold  on 

plot (states . x, abs (states . Zz) ) 
plot (states . x, real (states . Zz) , ’r ’ ) 
axis  tight 

borderlines (z , input) 

plot (states .x, abs (states . Zz) ) ;  plot (states .x, real (states . Zz) , ’ 
ylabel(’ Longitudinal  Stress’) 
subplot (222) 
hold  on 

plot (states . x, abs (states . Vz) ) 
plot (states . x, real (states . Vz) , ’r ’ ) 
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axis  tight 

borderlines (z , input) 

plot (states .x, abs (states . Vz) ) ;  plot (states . x, real (states . Vz) , ’r ’ ) 
ylabel(’ Longitudinal  Velocity’) 
subplot (223) 
hold  on 

plot (states . x, abs (states . Zx) ) 
plot (states . x, real (states . Zx) , ’r ’ ) 
axis  tight 

borderlines (z , input) 

plot (states . x, abs (states . Zx) ) ;  plot (states . x, real (states . Zx) , ’r ’ ) 
ylabel(’ Transverse  Stress’) 
subplot (224) 
hold  on 

plot (states . x, abs (states . Vx) ) 
plot (states .x, real (states .Vx) , ’r ’ ) 
axis  tight 

borderlines (z , input) 

plot (states . x, abs (states . Vx) ) ;  plot (states . x, real (states . Vx) , ’r ’ ) 
ylabel(’ Transverse  Velocity’) 

end 

end 

function  borderlines (z , input) 

[Ybounds]  =  get (gca, ’ YLim’ ) ; 

if  nargin<2 

for  j  =  l:length(z) 

plot([z(j)  z(j)] ,  Ybounds,’ — k’) 
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end 


else 

%make  rectangles  corresponding  to  types  instead! 
rheight  =  Ybounds(2)-Ybounds(l)+2; 
for  j  =  2 : (length(input . type) -1) 
rwidth  =  z(j)-z(j-l); 

if  input  .type (j )==’ s  ’  °/0orange  rectangle 

rectangle ( 'Position’ ,  [z ( j  —  1 ) , Ybounds ( 1 )  —  1 , rwidth, . . . 

rheight] , ’FaceColor ’, [255  177  86] /255, ’Linestyle’ , ’ — ’) 
elseif  input ,type(j )==’ f ’ 

rectangle ( ’Position’ ,  [z ( j  — 1 ) , Ybounds ( 1 )  —  1 , rwidth, . . . 

rheight] , ’FaceColor ’, [128  207  110] /255, ’Linestyle’ , ’ — ’) 

end 

end 

end 


end 
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APPENDIX  C:  MULTI-SOURCE  CODE 

The  code  below  is  the  code  for  the  model  presented  in  Chapter  4  for  generating 
plane  and  evanescent  wave  pressure  distributions  using  a  source  array.  The  model  is 
input  with  the  desired  frequency  and  angle  of  plane  wave  to  be  modeled  (an  evanescent 
wave  can  be  modeled  using  a  complex  angle  of  incidence),  along  with  the  length  of 
the  pressure  distribution,  the  standoff  of  the  array,  and  the  length  and  number  of 
sources  in  the  array.  It  outputs  a  vector  of  the  monopole  source  amplitudes  (which 
incorporate  phase  through  a  complex  component),  the  resulting  approximate  pressure 
and  velocity  generated,  the  correlation  coefficients  for  both  pressure  and  the  velocity, 
and  the  position,  pressure,  and  velocity  vectors  for  the  ideal  pressure  distribution. 

The  program  first  generates  the  ideal  plane  wave  pressure  and  velocity,  then  con¬ 
structs  the  matrix  for  calculation.  Once  the  complex  amplitudes  of  the  sources  has 
been  found,  the  program  calculates  the  pressure  and  velocity  resulting  from  the  mul¬ 
tiple  sources  and  the  correlation  coefficients  of  the  pressure  and  velocity  distributions. 

function  [A,Papprox,Vapprox,R_sq,X,Pideal,Videal,PdzApprox]  =. . . 
PV_approx (f , angle , L_x , ZO , L , N , phaseshif t ) 

if  nargin<7 

phaseshif t  =  0; 

end 

"^approximates  monopole  source  strength  using  least-squares  formulation 
%0UTPUTS : 

1  A  -  source  strengths 

°/„  Papprox  -  resulting  pressure  distribution  in  space 
°/„  Vapprox  -  resulting  velocity  distribution  in  space 
°/„  R_sq.P  -  fit  coefficient  for  the  pressure 
°/„  R_sq.V  -  fit  coefficient  for  the  velocity 
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%  X  -  x-coordinates  for 

*/.  Pideal  -  ideal  pressure  on  surface  (to  be  fit) 

*/.  Videal  -  ideal  velocity  on  surface 
7,  INPUTS : 

%  f  _  frequency 

7»  angle  -  complex  angle  of  incidence 

7»  L_x  -  length  of  desired  pressure  distribution 

7»  ZO  -  standoff  distance  of  array 

7»  L  -  length  of  array 

7»  N  -  number  of  sources 

7»  air  material  properties 
c  =  343; 
rho  =  1.2; 

7»  functions  of  frequency 
omega  =  2*pi*f ; 
k  =  omega/c; 

theta  =  real  (angle)  *pi/180+imag(angle)  *li ;  7«incident  angle  in  radians 

X  =  linspace  (-L_x/2  ,L_x/2 , 1024)  ’  ;  70position  vector 

7»ideal  pressure/velocity  distributions 

Pideal  =  exp(li*k*sin(theta)*X+phaseshift*pi/180) ’ ; 

Videal  =  -k*cos (theta) / (omega*rho)*exp(li*k*sin(theta)*X+phaseshift*pi/180) J ; 
7»  source  posiion  vector 

Y_x  =  linspace(-L/2,L/2,N) J  ;  7»  x-positions 
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Y  =  [Y_x  Z0*ones(N,  1)]  ;  “/.matrix  of  x,y  positions 

“/construct  green's  function  matrix 
G_P  =  zeros (length(X) ,N) ; 

for  i  =  l:length(X) 
for  j  =  1:N 

“/Green's  function  for  distance  r  from  source  to  point 
r  =  sqrt ( (X(i , 1) -Y( j , 1) ) ~2+(Y( j , 2) ) ~2) ; 

H(i,j)  =  exp(li*k*r)/(r) ; 


end 

end 

“/solve  system  HA  =  B  for  A 
B  =  Pideal ' ; 

A  =  pinv (H ' *H) *H ' *B ; 


“/„  generate  pressure,  velocity  from  sources 
Papprox  =  zeros (length(X) , 1) ; 

Vapprox  =  zeros (length(X) , 1) ; 

“/compute  approximate  pressure,  velocity  distributions 
for  j  =  1 :  N 

for  i  =  l:length(X) 


r (i)  =  sqrt ( (X(i , 1)-Y(j , 1) ) ~2+(Y(j ,2) ) ~2) ; 
phi (i)  =  atan( (X(i , 1) -Y(j , 1) )/Y( j , 2) ) ; 
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Papprox(i)  =  Papprox(i)+A(j)*exp(li*k*r(i))/(r(i)) ; 
Vapprox(i)  =  Vapprox(i)-A(j)*l/(rho*c)* . . . 

(l-l/(li*k*r (i) ) ) *exp(li*k*r (i) )/(r(i))*cos (phi (i) ) ; 


end 

end 

°/0calculate  R  from  P,V 

MSE. P=  sqrt (sum( (real (Pideal)-real (Papprox’ ) ) . ~2) ) ; 
RMS.P  =  sqrt (sum(real (Pideal) . “2) ) ; 

MSE.V  =  sqrt (sum ( (real (Videal) -real (Vapprox’ )). ~2) ) ; 
RMS.V  =  sqrt (sum(real (Videal) . "2) ) ; 


R_sq.P  =  1-MSE. P/RMS. P; 
R_sq. V  =  1-MSE. V/RMS .V; 


VITA 
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